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Abstract 


Chemical vapour deposition (CVD) reactors are used for growing thin films of 
special materials on suitable substrates. A particular reactor at Centre for Ad- 
vanced Technology, Indore is in the development stages to manufacture dome- 
shaped infrared sensors for night vision applications in the defence industry. The 
computational study of transport phenomena in such a reactor has been accom- 
plished in the present work. The reactor has a central jet carrying hydrogen 
sulphide and peripheral jets carrying zinc that react in the core to form zinc 
sulphide as a deposit on the solid surfaces. The reactants are transported by 
argon that acts as an inert carrier gas. The primary deposition surface is a trans- 
versely placed block within the reactor tube. A two-dimensional model has been 
developed to study flow and mass transfer phenomena within the reactor. Three 
geometric configurations of the reactor have been studied and two new geometries 
have been proposed with the aim of providing alternatives to the proposed CVD 
technology. 

The mathematical model comprises of the system of conservation of mass, 
momentum and the species transport equations. The spatial temperature varia- 
tion in the reaction zone is expected to be small and has been neglected. The gas 
phase reaction of zinc and hydrogen sulphide is taken to occur in a single step, 
besides being homogeneous and irreversible. The dilute solution approximation 
has been uniformly adopted in the present work. Higher order phenomena such 
- as the thermo-diffusion and dufour effect have not been considered. 

The momentum equations have been solved by finite volume method (FVM). 
The dependent variables have been arranged on a collocated grid. Momentum 
interpolation has been used to resolve the difficulty of decoupling velocity and 
pressure. For the average pressure and temperature encountered in the reactor, 
Reynolds numbers are in the range of 10-100. Hence, laminar incompressible con- 
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ditions are realized in the flow field. The mass transfer equations have been solved 
by the application of the operator splitting algorithm. The source term that mod- 
els chemical reactions in the species transport equation has been linearized with 
respect to the species concentration. This permits a partial analytical solution 
to the governing equations and improves the computational efficiency. The de- 
position surfaces have been modelled to be passive. It was seen that arriving 
at a correct boundary conditions for a passive surface was difficult and led to 
errors in the deposition rate profiles. This situation has mitigated by employing 
a Dirichlet condition on the substrate that slightly overestimates the deposition 
rate. 

The following results have been obtained in the present work. The species 
concentration field correlates closely with the flow field in the reactor. Over the 
range of Reynolds numbers studied, the flow pattern reveals a recirculation zone 
adjacent to the side walls and a stagnation zone ahead of the block. Species mix- 
ing predominantly occurs in the zone closer to the reactor axis and the products 
are carried away by the stream bypassing the block. Under favourable condi- 
tions, a considerable amount of deposition can occur over the block surface. The 
performance of the reactor designs proposed at Centre for Advanced Technology, 
Indore are not satisfactory in terms of deposition rate and its uniformity. The 
two new geometries explored in the present work not only improve the deposition 
rate by a factor of two (1.027 jim/min from 0.542 fim/min) but also retain the 
uniformity in the thickness of the layer formed. The new geometries have the 
potential to significantly reduce overall cost and time (from over a year to just a 
few months) for thin film deposition. 
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Chapter 1 


Introduction 


Thin films have been the topic of a large number of investigations during 
the last two decades. Thin films are technologically important, particularly in 
the fields of semiconductor electronics and computer hardware. Thin films can be 
prepared by using a variety of methods, among which chemical vapour deposition 
is the most evolved and has received widespread acceptance. Chemical Vapour 
Deposition (CVD) refers to the technique of growing an epitaxial layer on a sub- 
strate by reactions in chemical species in the vapour phase across an activation 
energy barrier. 

Applications of the CVD technique can be seen in different engineering sys- 
tems, most importantly in microelectronics fabrication. It is used to produce 
highly uniform thin layers (0.01 — 10 /j,xn) of semiconductors such as epitaxial 
silicon and gallium arsenide, dielectrics such as silicon dioxide and silicon nitride. 
Other applications include the protective coatings of AI2O3, TiC, SiC and B4C; 
anticorrosive coatings such as BiN, MoSi2, SiC and B4C; fibers such as B, B4C 
and SiC for composite materials and high-purity monolithic compounds such as 
ZnSe, ZnS, CdS and CdTe as infrared sensors. 

The basic components of a CVD reactor are shown in Figure 1.1. It mainly 
contains nozzles in the inlet plane for the injection of reactants, a curved sub- 
strate, a zinc dump box where the exhaust is collected. The processes inside 
a CVD reactor follow a sequence. These are outlined with the help of Figure 
1.2. The incoming reactants are carried to the reaction zone by the basic trans- 
port mechanisms occurring inside the CVD reactor. Here they undergo a series 
of chemical reactions that generate product species. The newly formed species 
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main gas flow 


exhaust of the reaction products 
fiom the reactor 



adsorption of surface surface reactions, 

film precursor diffusion nucleation and 

lattice incorporation 


Figure 1.2: Schematic representation of the sequential processes occur in CVD 
reactors. 

along with the remaining reactants are then transported to the susceptor. The 
physical species present at this location are adsorbed to the surface and are later 
mobilized by surface diffusion. The reactant species can also react with the sus- 
ceptor and intermediates can be formed. The intermediates are either adsorbed 
or react to form a more stable product. In addition, nucleation and lattice incor- 
poration lead to a thin layer of product species on the susceptor. The gaseous 
products desorb at the solid surface and diffuse away from the susceptor to the 
main gas flow. Lastly the gaseous products are exhausted with the main gas flow 
from the reactor. 

CVD reactors may operate at atmospheric pressures (APCVD), or at low 
pressures (LPCVD). The operating pressure for the first type generally varies be- 
tween one-tenth of an atmosphere to an atmosphere (0.1 — 1 atm). In the latter, 
pressures are typically in the range of (10"^ — 10“^) atm. CVD reactors are also 
classified by the source from where they derive the energy for the reactions to 
proceed. The thermally activated reactor has gained importance over the years 
due to two factors. First, it provides versatility and secondly a high through- 
put. In this category plasma-enhanced CVD (PECVD), photo-CVD(PCVD) and 
electron-beam assisted CVD are the most important. Clearly the steps men- 
tioned in Figure 1.2 are temperature dependent. At lower temperatures, surface 
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reactions are not fast and they constitute the rate limiting step. At higher tem- 
peratures surface reactions are quite fast compared to the mass diifusion from 
and to the surface. Thus at the higher temperature it is diffusion between the 
solid and gaseous phases that limits the transport mechanism. 

Over the years various CVD geometries have evolved, horizontal reactor, 
pan-cake reactor and barrel reactor being the most explored. An excellent com- 
pilation of industrial reactor geometries have been given by Mahajan (1996). 
With the continuing trend towards larger-diameter wafers, higher quality and a 
decrease in minimum feature size in microelectronics manufacturing, CVD tech- 
niques will call for better understanding of the underlying transport phenomena 
and its relationship with the deposition characteristics. It will be a challenge to 
achieve ±2% unifomity in layer thickness of 0.18 nm. From Figure 1.2, it is clear 
that the underlying transport phenomena during chemical vapour deposition is 
quite important and has a strong influence in the deposition layer thickness, its 
uniformity and purity. 


1.1 Literature Survey 


A number of investigations comprising numerical simulations and experi- 
ments have been carried out over the last decade. Most of these investigate the 
role of transport phenomena in the overall growth characteristics in a CVD re- 
actor. Growth of epitaxial silicon layer has been focused on by most authors. 
Emphasis has been given to the modeling aspects of the transport mechanisms 
inside the reactor. Selected papers on numerical modeling of CVD reactors are 
reviewed in the present section. 

Moffat and Jensen (1988) studied the effect of three dimensional flow phe- 
nomena on critical parameters such as the deposition rate on the substrate sur- 
face. Species distribution in the reactor and the onset of natural convection in a 
horizontal reactor have also been studied. A system of three dimensional equa- 
tions was solved numerically in an axisymmetric domain with the reactive mass 
transfer problem being decoupled from momentum and energy, by the assump- 
tion of a dilute solution. Applicability of a similarity solution and its limitations 
have been pointed out. Homogeneous chemical kinetics with surface reactions 
have been considered. The mechanical pressure in the reactor has been split in 
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the form 

P{x, y, z) = P{x) + v\x, y, z) 


where P denotes the bulk pressure used in the axial momentum equation, driving 
the strong axial convective flow. The perturbation p' denotes the variation about 
the mean pressure and drives the transverse flow. It has been emphasized that 
temperature boundary conditions on the side wall influence the distribution of 
the species. Two types of boundary conditions have been used: (1) adiabatic, 
which is applicable for an air cooled system and (2) side wall temperature equal 
to that of the top wall, relevant in water cooled systems. In the base line study, 
1 atm has been taken as the mean reactor pressure. 

Klejin et al. (1989) solved the flow, energy and the species transport equa- 
tions in two dimensional form in an axisymmetric geometry. The authors strongly 
recommended that the coupled system of equations should be adopted to analyze 
the complex transport phenomena inside a CVD reactor. Thermo-diffusion effect 
and multicomponent diffusion were taken into account. One interesting aspect 
revealed was that a deviation from the impinging jet arrangement to stagnation 
flow improves the growth uniformity. It was pointed out that to keep the recircu- 
lation flow caused by natural convection suppressed higher flow rates and lower 
total pressure are favorable. Even at 10 torr pressure, recirculation patterns ap- 
pear owing to free convection. 

Takoudis et al. (1991) developed a mathematical model for epitaxial silicon 
growth in a pancake reactor. The equations governing momentum, energy and 
species transport were numerically solved. Free jet flow was assumed due to the 
jet-to-reactor diameter ratio being smaller than 5. Both homogeneous gas phase 
reaction and the complex surface chemistries were modeled. It was pointed out 
that if molecular weights of the reactants and product species are quite differ- 
ent, thermo-diffusion effect would be important. This term contains the product 
of the species concentration. In the dilute solution limit it becomes low. The 
authors showed that with increase in flow rates and decrease in the susceptor 
Damkohler number, the deposition layer uniformity improve. 

Duverneuil and Couderc (1992) presented a two dimensional model for sil- 
icon LPCVD. A sub-atmospheric pressure of 0.1 torr was employed for the re- 
actor. Governing equations have been solved by the stream function- vorticity 
method taking a representative inter-wafer region as the physical domain. It was 
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hypothesized that a large number of phenomena linked to the flow, transport pro- 
cesses, chemical reaction reach dynamic equilibrium. They reproduce an identical 
scenario from one inter-wafer spacing to the next. Boundary conditions of the 
periodicity type were explored and analyzed in view of the ease of computations. 
Two-step gas phase reaction of silane pyrolysis and silylene insertion were taken 
as the chemical reaction model. By comparing the results of the study with a full 
two dimensional gas phase reaction model it was shown that in the latter, growth 
rate is overestimated due to the absence of rapid surface reactions. 


1.2 Statement of the problem 



In the present study the reactor that has been analyzed does not match 
the CVD geometries reported in the literature. It has the characteristics of both 
the stagnation flow and the impinging-jet LPCVD Mahajan (1996) reactor. A 
vertical cross section and the side view of the reactor are shown in Figure 1.3. 
The basic geometry consists of a central jet and four peripheral jets positioned at 
an angle of 90° with respect to each other. Argon (A), the carrier gas mixed with 
hydrogen sulphide (H 2 S) is injected through the central jet while argon carrying 
zinc vapour (Zn) is introduced through the peripheral jets. For the simulation a 
two dimensional plane consisting of the central jet and two peripheral jets such 
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Figure 1.4: The axisymmetric two-dimensional reactor geometry employed in the 
present study with a fiat (a) and a curved (b) deposition surface 




Figure 1.5: The axisymmetric two-dimensional reactor geometry employed in the 
present study with a convex deposition surface and jets separated by a distance 
(a) concave deposition surface and coaxial jets (b) 
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as plane AA' of Figure 1.3 is considered. Due to the inherent symmetry of the 
reactor, an axisymmetric configuration has been considered in the present study. 
Two deposition surfaces have been considered as in Figure 1.4. Figure 1.4(a) 
shows the deposition surface to be fiat and Figure 1.4(b) shows the deposition 
surface to be a hemispherical bowl. The substrate surface is concave for the 
geometry shown in Figure 1.4(b). A variant of this geometry where the surface 
is convex in shape (1.5a) has also been addressed briefly. With the surface as 
concave but the jets being coaxial (1.5b) has also been studied. 

The CVD reactor described above has been proposed to be built at CAT, 
Indore. The operating parameters for the proposed reactor are as follows; 


1. Mass flow rate of argon through the central jet = 500-600 cc/min 

2. Mass flow rate of argon through the peripheral jet = 37.5-50 cc/min 

3. Vaporization rate of zinc through the peripheral jet = 7-10 gm/hr 

(mass fraction of zinc at the inlet is approximately 0.08 which is the ratio of mass 
flow rate of zinc through the peripheral jet and the net mass flow rate of argon) 

4. Mass flow rate of hydrogen sulphide = 90-100 cc/min 

(mass fraction of hydrogen sulphide at the inlet is approximately 0.1 which is the 
ratio of mass flow rate of hydrogen sulphide through the central jet and the net 
mass flow rate of argon) 

5. Mean reactor pressure = 30 mbar 

6. Substrate temperature = 900 K (maximum temperature 1023 K) 

7. Temperature at the central zone of the reactor = 923 K 

8. Temperature of the outer most insulation at the reactor wall = 323 K 


1.3 Thesis Organization 


In Chapter 2 the set of governing equations, boundary conditions with all 
the approximations are written and discussed. Several modeling aspects relevant 
in the present study are discussed there which help to formulate these equations. 
Numerical techniques used to solve the governing equations and applicable bound- 
ary conditions are discussed in Chapter 3. Solution algorithm for momentum and 
mass transfer have been outlined. Results of grid independence test have been 
presented in Chapter 3. Results for the flat and curved substrate geometries have 
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been presented in Chapter 4 followed by results of the two alternative geometries 
in the same Chapter. In appendix A calculations for the species properties have 
been shown briefly. Validation of the fluid flow and the mass transfer codes have 
been shown in appendix B. The derivation of orthogonal mapping equations used 
in grid generation have been outlined in appendix C. 

1.4 Objectives of the present work 


The main objective of the present work has been to develop a mathematical 
model to study the various aspects of transport mechanism inside the proposed 
CVD reactor. The other objectives are listed as below: 

1. To build a clear understanding of the implications of the various modeling 
approximations on the overall deposition process. 

2. To report the computational artifacts related with the mathematical modeling. 
Specifically the modeling differences between an chemically active and a passive 
substrate. 

3. Achieve a steady state solutions for the flow and concentration fields arising 
during the growth process. Analyze these solutions in the view of deposition 
characteristics. 
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Chapter 2 


Mathematical Formulation 


The chemical vapour deposition process is intricately controlled by fluid flow 
and transport phenomena in the reactor. The equations that govern the velocity, 
temperature and species concentration fields enforce the physical principles of 
conservation of mass, Newton’s second law of motion and conservation of energy. 
These must be supplemented by equations that model chemical reactions among 
the species present in the carrier gas. The assumptions that lead to the reduced 
set of equations applicable for the proposed CVD reactor have been discussed. 
Several modeling aspects required in this connection have been discussed in detail. 
The mathematical formulation presented below draws heavily from the earlier 
works of Moffat and Jensen (1988), Klejin et al. (1989), Takoudis et al. (1991), 
Duverneuil and Couderc (2000). 


2.1 Equations Governing Flow, Heat and Mass 
Transfer 

2.1.1 Continuity equation 

The continuity equation expresses the conservation of mass in the fluid flow do- 
main. It is written below in coordinate free form : 

^ + V.(H = 0 (2.1) 
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Expanding and writing it in a cartesian coordinate system, we get 


^ djpu) d{pv) djpw) _ 

dt dx dy dz 


( 2 . 2 ) 


In the present study, the fluid flow inside the CVD reactor is limited to low ve- 
locities. Moreover, the carrier gas which carries the reactant gases to the reaction 
zone and the product gas to the deposition zone is in great excess, compared 
to the other species, leading to a dilute solution. Thus density variation inside 
the reactor can be neglected and the fluid flow can be considered to occur in the 
laminar, incompressible regime. Hence equation (2.1) reduces to 


V.u = 0 


which when expanded takes the form 

du dv dvj 
dx'^ dy'^ dz 

2.1.2 Momentum Equations 


Navier-Stokes equations govern the transport of momentum in the fluid flow do- 
main. For incompressible, laminar flow, they can be written in compact form 
as 

,au „ , „ 


p(^ -I- u.Vu) = -Vp -h pV^u 


Here the density and the viscosity of the carrier gas have been taken to be constant 
and body force to be negligible. It is convenient to combine the left side of 
equation (2.3) with the incompressible condition V.u = 0 to yield the so called 
conservative form: 

,d\x „ , ,au „ , 

+ u.Vu) = p(— + V.(uu)) 


The three momentum equations in cartesian coordinates are written in conserva- 
tive form as follows : 
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x-momentum eqn. 


.du d{v?) d{uv) d{uw). 


dp .d'^u d'^u d'^u 

al ^ ^ V ^ 


(2.4) 


y-momentum eqn 

.dv divu) d{v^) d(vw). 

p(i57 + 4r^ + -^ + -V) 


' dt dx 


dp d'^v d^v d'^v 
dy ^ dx'^ dy'^ dz^ 


(2.5) 


z-momentum eqn. 

.dw d{wu) d{wv) d{w‘^). 


dp d‘^w d^w d'^w 
dz ^ dx"^ 9y^ dz‘^ 


( 2 . 6 ) 


2.1.3 Energy Equation 


The equation that governs the transport of energy in a fluid system can be ob- 
tained form the conservation of energ}'. In cartesian coordinates, it can be written 
as 

^ ^ ^ ^ ^ 

dt dx dy dz dx'^ dy‘^ dz^ 


The following assumptions have been made in deriving the energy equation in 
the form given by equation (2.7) 


1. Thermal diffusivity does not vary considerably with the temperature in the 
sense that an average value can be used in the analysis. This is acceptable because 
the CVD reactor operates under nearly isothermal conditions. 


2. There is no heat sources in the flow domain. Though the reaction taking place 
inside the reactor is exothermic, the heat release has been taken to be minor. 


3. Carrier gas is in considerable excess compared to the reacting species leading 
to a dilute solution. Thus concentration gradients do not lead to heat fluxes in 
the field. In other word, the dufour effect, also called the diffusion-thermo effect 
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has been neglected in the present study. It has been shown that this effect is 
only of second order importance in practical CVD systems Bird et al. (1960) and 
Hirschfelder et al. (1964). 


4. As velocities involved in the transport processes are not high and consequently 
viscous dissipation effect can be neglected. It is well established that the dissipa- 
tive effect due to the internal shear is important only in cases where the velocities 
are super-sonic or nearly so. 


2.1.4 Species transport Equation 

In thermally activated CVD reactors, we frequently encounter species transport 
processes coupled with chemical reactions. The governing equations for species 
transport is written here for the ith species using diffusive mass fluxes relative 
to mass-averaged velocity of the gas mixture. The mass-averaged velocity can 
be obtained by solving the Navier-Stokes equation and the resulting equation is 
similar in form to the energy equation and to each component of the Navier-Stokes 
equations. Let oij be the mass fraction of ith species. The transport equation for 
each of the i species can be written as Mahajan (1996) 

= + (2.8) 

i=i 


where the total diffusive mass flux Jj of species i is composed of Jf, the flux due 
to ordinary diffusion caused by a concentration gradient, the pressure diffusion 
Jf , the forced diffusion jf and the thermal diffusion Jj. Hence 

ji = j; + jf + jf + jr (2.9) 


The ordinary diffusive flux in a multicomponent mixture is given by the Stefan- 
Maxwell equations Bird et al (1960). In terms of mass fractions and species 
fluxes, they are written for a system containing n components as 

•m ^ Tl ^ 

Vui+w.V(liiM) = - 
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where M is the average molecular weight of the mixture given by 

n 

M = (2.11) 

i=i 


Here Xj denotes mole fraction of the 2 th species. Equation (2.10) is coupled to 
the steady state mass balance constraint 

n 

= “ P-12) 

2=1 


The above equations 2.10 - 2.12 form a closed set of equations for the mass fraction 
Wj. The n diffusive fluxes can then be obtained as 



n 

-pDeiVu)^ - poJ^DeiV {In M) + MuJ^Dei 



(2.13) 


where Dei the effective diffusion coefficient for species i is given by 


Dei = 


n 

(E 


MjDi, ^ 


(2.14) 


Equations (2.10) and (2.12) can be solved iteratively in conjunction with equa- 
tions (2.13) and (2.14), to obtain the ordinary diffusive flux for all the species. 

The thermo-diffusion term Jf , which accounts for the diffusion of a species 
due to temperature gradient, is given by 

Jf = -DfV(lnT) (2.15) 


where Df is the thermo-diffusion coefficient for the species i. An approximate 
expression for the thermo-diffusion coefficient is given by Bird et al. (1960), 
Hirschfelder et al. (1964) 

Dr= E -MiMiDijkl ( 2 . 16 ) 
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The quantity is the thermo-diffusion ratio which varies from one mixture to 
the next. In a majority of applications, the lighter species drifts towards the zone 
of high temperatures and the heavier species drifts towards the zone of low tem- 
perature. This phenomena is called the Soret effect. It can be seen from equations 
(2.15) and (2.16) that this term is coupled with the species concentration. For a 
dilute solution the individual species concentration levels are quite low and this 
term is only of secondary importance. This limiting condition has been enforced 
in the present work. Details of the thermo-diffusion effect have been discussed 
by Hirschfelder et al (1964). 

The last term in the species transport equation (2.8) expresses 

the creation or consumption of species due to a combination of gas phase chemi- 
cal reactions. The rate expressions for the chemical reactions involved depend on 
reaction kinetics and stoichiometry. 

The present work is restricted to modeling elementary, homogeneous reac- 
tions in CVD reactors. Elementary reactions proceed in a single step and their 
rate expressions directly employ the stoichiometric coefficients. It is difficult to 
obtain the rate expression for nonelementary reactions because of the appearance 
of intermediate molecules. Nonelementary reactions account for the fact that a 
single reaction is in reality the overall effect of a sequence of several elementary 
reactions. In the present work, the intermediates are assumed to be so little 
in amount that they escape detection. The second useful scheme for classify- 
ing chemical reactions is to identify them as homogeneous or heterogeneous. A 
homogeneous reaction model involves only one phase of the species, usually the 
gas phase. For a heterogeneous reaction the species are required to be present 
in at least two or more phases. These two types of reactions have significant 
implications in mathematical modeling. A homogeneous model accounts for pro- 
duction/destruction rates of species throughout the domain while a heterogeneous 
reaction is restricted to a small portion of the flow domain. For CVD systems 
that produce silicon substrates this zone is the susceptor boundary. Much of the 
research prior on CVD systems has employed a heterogeneous reaction model 
with a focus on surface reactions. 

In the present study, a CVD reactor in which elementary (single-step), 
homogeneous reaction takes place in the gas phase has been modeled. The sub- 
strate surface mounted normal to the main flow direction is chemically inert. The 
deposition is neither due to surface reactions nor is surface diffusion important. 
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Moreover the reaction studied is instantaneous, in the sense of large reaction con- 
stants, and irreversible. Thus all reactions are confined to the gas phase. The 
issue of deposition has been resolved by examining the physical nature of the reac- 
tion and its products. The product species is expected to be composed of micron 
level particles and so the surface forces exerted on these particles overwhelm their 
weight. Thus the product species is assumed to be transported by the carrier gas 
to the substrate where it is deposited at a rate fixed by the fiuid-to-solid diffu- 
sion fluxes. Thus the particles formed away from the deposition surface can be 
considered to be in gas phase driven by the base flow field towards the deposition 
surface. On the deposition surface the product species can be considered as a 
solid while the reaction takes place in the gas phase. In this respect, considering 
the reaction to be homogeneous is justified. As the product species deposited 
on the substrate does not return to the flow, the reverse reaction has not been 
accounted for. 

In view of the above approximations, the stoichiometry involved in the 
present work is quite simple and a single step elementary reaction has been cho- 
sen to express the rate kinetics. The rate expression of a chemical reaction acts 
as a mass source of species in the transport equation, but depends on the fields of 
concentration and temperature. The temperature dependence is due to the fact 
that reaction rate constants obey the Arrhenius’ law Levenspiel (1999) given by 

k = A:oexp(— -^) 

where ko is the pre-exponential factor, Ea is the activation energy and R is the 
universal gas constant. In the present study an arbitrary large value compatible 
with double precision arithmetic has been assigned to the rate constant as the 
reaction is assumed to be instantaneous in nature. 

Consider a reversible reaction of the type 

UlTi + U2r2 H KiPi -I- K2P2 -I 

k-r 

where r and p denote reactants and products, respectively. The net rate of 
creation of species j is given by 
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where and 3?^ are the forward and backward reaction rates, respectively, given 
by 


I 

= kfY[[P^n) 

I 

= ^rj][cp, 

I 

= hWipf^v^) 


Thus for a reaction containing four species 


vaA + vbB i/qC + i^dD 

kr 

(2.17) 

the reaction rates for species C is given by 


3?c = 

(2.18) 

where the forward and backward reaction rates are given by 


= kf{pUJA){pOJB) 

(2.19) 

= kr{pUJc){pUJD) 

(2.20) 

For an elementary (single-step) reaction, the rates of reactants and products are 
related by the stoichiometric coefficients in the following relationship Levenspiel 

('1999') 

(2.21) 



Several additional simplifications appropriate for the CVD reactor being modeled 
are listed below 
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1. In atmospheric and sub-atmospheric CVD reactors, the carrier gas is in con- 
siderable excess compared to the reactant species in order to avoid homogeneous 
nucleation of particulates. Thus the dilute solution approximation can be en- 
forced to develop the transport equation. 


2. In view of the dilute solution approximation, the presence of other species on a 
certain transported species can be neglected. Thus ordinary diffusion is governed 
by the Pick’s law of binary diffusion. In case of a multicomponent diffusion system 
equations (2.10) and (2.12) are solved to get the ordinary diffusive fluxes of each 
species. These are then inserted in the transport equation to solve for the species 
concentration field. 


3. In a low pressure CVD reactor, the pressure gradient encountered is not 
significant enough to separate the chemical species. Hence the pressure diffusion 
term can be neglected. 


4. The forced diffusion term is identically zero if gravity is the only external force 
Mahajan (1996). 


5. The thermo-diffusion effect may be considerable in regions of sharp tempera- 
ture gradients and in multicomponent mixtures where the molecular weights of 
the species are not close to one another. For the C\'D reactor under discussion, 
temperature changes about the mean value of 700° C are within ±1° C. Hence 
the thermo-diffusion effect can be neglected in the formulation. 


6. In view of the near constancy of temperature in the reactor, temperature 
dependence of the binary diffusion coefficients and density of the carrier gas have 
been neglected. 


With the help of above mentioned approximations, the diffusive mass flux for the 
ith species is given by 

3i = ri = -DijV{poj,) ( 2 . 22 ) 
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with j corresponds to the carrier gas and Jf is the ordinary diffusion flux of 
species i into the carrier gas. Hence the simplified form of the species transport 
equation takes the form 

+ V.(pua;0 = V.(- A; V(pa;0) + ^ (2.23) 

By inserting the rate expressions, namely equations (2.19) and (2.20), into the 
governing equation we can form the complete set of equations that govern the 
species concentration field. For a reaction of the form given by equation (2.17), 
the governing equations have been written explicitly for all the species as follows 


species A 

d{pUJA) , d{pUUlA) , d{pVUlA) , d{pWUJA) _ p, ,d‘^U)A , d'^UJA , d'^COA. 

dt dx dy dz ' ^ dx^ dy^ ^ dz'^ ’ 

+UAMA[kr{puJc){p0JD) “ A:/(pw^) (pws)] (2.24) 

species B 

d{pUB) djpUOJB) djpVOJB) djpWUJB) , d'^U)B d'^OJB d'^OJB . 

dt dx dy dz P b[ Qy2 

+UBMB[kr{pUJc){p(BD) “ kf{pU)A){pUJB)] (2.25) 


species C 


d{pujc) , djpucxc) , d{pVLOc) , djpwojc) _ ^ ^ d'^ug , d^uc , d^ujg ^ 

dt dx dy dz ^ ^ dx"^ dy'^ dz^ 


+i/cMc[A:/(pwa)(pwb) - ^r(pwc)(pw£))] (2.26) 


species D 

d{pOJD) 

dt 


d^pnojo) d{pvu)D) d^pwixo) 
H X \ X b 


dx ' dy ' dz dx'^ 

+UDMD[kf{pU}A){p(jOB) - kr{p(Xc){pUJD)] 


fd'^ujr, ^ d^ojD , d'^ujD-, 
P^D[ o_ 2 ”b 0..2 Xl~2 ) 




dz"^ 
(2.27) 
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where Ui corresponds to the stoichiometric coefficient of the zth species in the 
chemical reaction given by equation (2.17). 


2.2 Nondimensionalization of the governing equa- 
tions 


Nondimensionalization of the governing equations facilitates parametric study 
over a smaller set of flow properties and geometric variables to get an insight into 
the transport phenomena inside the CVD reactor. The scales for the dependent 
and independent variables for normalization are written below : 


u = 


u 

* 'y 

w 


= -£- p*-. 

p 


tUa 

* ^ 

Uc' 

3 

II 


> P 

Pea 

~ pur 

t = 

d ’ 


y 

* ^ ^ ^ 

T- 

Tb 

, Wi 

CJ, = , 

SR! = 


, D* 

_ 

d’ 

d’ 

Ta- 

Tb’ 

% ' 

3 

P^Xo 

7 Z 

Di,o 


Here Ta and Tj are two convenient reference temperatures in the reactor. Sub- 
stituting the normalized quantities in the governing equations we get the set 
of nondimensionalized equations. These are presented below after dropping the 
superscript (*). 


continuity equation 


dp d{pu) d{pv) d{pw) 

dt dx dy dz 


(2.28) 


X momentum equation 

dipu) div?) d(uv) diuw) dp 1 .d'^u d'^u d'^u, ... 


y momentum equation 


d{pv) 


+ 


d{vu) d{v^) d{vw) 

dx dy dz 


dp 

/ + 

dy 


1 d'^v 
Re dx'^ 


d'^v 




dt 


(2.30) 



z momentum equation 


d(pw) 

dt 


+ 


d{wu) 

dx 


d{'wv) 

dy 


+ 


d{w^) 

dz 


dp 1 .d'^w d'^vj d‘^w. . , 


where Re the Reynolds number is based on the central jet velocity, 

reactor diameter and the viscosity of the carrier gas at the average reactor tem- 
perature. 


Energy equation 

^ am ^ ^ 

dt dy dy dz Pej^ dx'^ dy'^ dz'^ 

T I rl 

where Pej^ = the Peclet number for heat transfer is based on the central jet 
velocity, diameter of the reactor and the thermal diffusivity of the carrier gas at 
the average reactor temperature. 


Species transport equations 


species A 


d{pU)A) , d{pUijJA) , d{pV(jJA) , d{pWUJA) 
H 1 r 


dt 


dx 


+ 


dy 


dy 

{pDa- 


dz 


1 Da,o ( d ^ ^ dojA 
FemDA/dx^^ ^ dx '' 


d , _ dujA. d , ^ dujA 


dy ^ ^ dz^^^^ dz 


)] 


d Pc Ma,, , ,, 

+^/l ~^A,o {kr<^C<^D - kfOJA(^B) 


■Ur. 


(2.33) 


species B 


d{pUJB) , d{pUU)B) , d{pVU}B) , d{pWUB) _ 1 DB,o;d duJB-^ 


dt 


dx 


dy 


dz Pern 0^5^ 
dujB^ 


d / r-\ ^ / T~\ M 

+Vb Jj~^A, 0 ^ jJ- [kr^C^D - kfOJAOJB) 


(2.34) 
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species C 
djpcoc) 
dt 


+ 


d{puuc) d{pvuic) d{pwujc) 

H r h 


1 D, 


C.o 


dx dy dz 

d dujc. 


Pem D 


m 


d dujc.. 


+ - krUciXlD) 


.5 duJc. 


(2.35) 


species D 


d{piOD) d{puuD) dipvuo) d{pwujD) 

l” r\ i-v 


dt 


dx 


dy 


dz 

d 


1 Defied diXD 

VtraDA^dx^^ ^ dx ^ 


,d duD. d du)D.. 

^dy^^ ° dy^^ dz^^ ^ dz 

+ l^D^^A,o-^^{kfU}AOJB - KuJcOJd) 


Mr Mr 


(2.36) 


where Pem = irf^ Peclet number associated with the mass transfer based 
on central jet velocity, diameter of the reactor and the binary diffusion coefficient 
of species A into the carrier gas at the average reactor temperature. 

In the present study, the system of equations including the continuity equa- 
tion, momentum equations and the species transport equations are taken to be 
segregated in the sense that the flow field is independent of the concentration and 
temperature fields. This is justified because the reactants and products are in a 
dilute solution. Moreover, with the coupled system approach, the computational 
effort is enormous. Any effect pertaining to the momentum transport and the 
species transport related to the changes in the temperature field have not been 
studied. Hence the energy equation has been excluded from the governing equa- 
tions. The focus of the present work is on the steady state concentration fields 
of species A - D and the deposition rate on the substrate placed in front of the 
incoming jets. 


2.3 Initial and boundary conditions 


The governing equations of section 2.2 require both initial and boundary 
conditions for a complete solution. These are summarized below. 
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2.3.1 Initial condition 


The initial condition is quite important from a computational point-of-view. How- 
ever when a steady state solution is sought, any initial field can be assigned to 
the dependent variables. In addition, converged fields at a lower Peclet number 
or a coarser grid can also be employed. These aspects are discussed m chapter 3. 


2.3.2 Boundary conditions 

Boundary conditions at various surfaces of the reactor are discussed here 


Solid surfaces 


At solid surfaces, no-slip conditions are applicable for the three components of 
velocities, written as 

14 = 0, u = 0, w = 0 (2.37) 


For the energy equation either constant wall temperature or constant wall heat 
flux are prescribed at a solid surface, namely 

T = specified or q = — fcVT = specified (2.38) 


For the species transport equations, the net mass flux of fth species must equal 
the rate of production of species i through rij surface reactions at the surface. 
Hence 

ris 

(pcjiu -h Ji).n = (2.39) 

j=i 


With the absence of any surface reactions and the no-slip velocity boundary 
conditions the above condition reduces to the following form with only ordinary 
diffusion in place : 

Ji.n = J^n = 0 (2.40) 
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where n is the normal unit vector on the surface. In the present study, all the 
solid surfaces are passive and impermeable. Hence the above boundary condition 
(2.40) is appropriate for the species transport equations. This condition have 
special significance for the deposition surface and is discussed in detail in chapter 
3. 

Inlet plane 

At the inlet, velocities are prescribed as 

u, V, w = specified (2 41) 

for the energy equation temperature is specified at the inlet ; 

T = given (2.42) 

For the species transport equation prescribed species concentration and no net 
diffusion flux at the walls are appropriate. The species concentration can be 
calculated as the ratio of known mass flow rates of the various species and the 
carrier gas. Hence, it follows that 

Wj = specified at jet locations and .n = 0 on solid surfaces (2.43) 

Exit plane 

At the exit plane sufficiently away from the reaction and deposition zone, zero 
derivatives of the dependent variables in the stream-wise direction has been en- 
forced : 

V(f).n = 0 (2.44) 

where n is the unit vector along the stream-wise direction at the exit plane. Here 
(f) denotes all the dependent variables such as velocities u,v and w, temperature 
T and concentration Wj. 
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At the Axis 


For the reactors with axisymmetric geometry no radial variation of the dependent 
variables at the axis is permitted. Hence 

(V<?i').n = 0 (2.45) 


where (j) denotes all the dependent variables namely velocities, temperature and 
concentration and n denotes the unit normal vector on the axis. 


2.4 Range of parameters 


To get a clear idea about the transport mechanism inside the reactor and the 
effect of geometry, a parametric study has been conducted. While carrying out 
the study one parameter has been changed at a time, keeping all other variables 
fixed. The effect of this parameter on the transport phenomena and the deposition 
rate on the substrate inside the reactor have been studied. The dimensionless 
parameters relevant to the present work are written below : 


Reynolds number 


In fluid mechanics Reynolds number signifies the ratio of the inertia force and 
the viscous force acting in tandem, given by 

where Uref is the reference velocity of the problem considered, Lc is the character- 
istic length associated with the fluid flow, p and p, are the density and viscosity of 
the fluid. At various Reynolds numbers we get an idea about the viscous effects 
on the flow and the size of the recirculation region in relation to the position of 
the substrate. It also gives us an idea of the effect of mass flow rate on the depo- 
sition profile because Reynolds number is directly linked with the mass flow rate 
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of the carrier gas. In the present study the range of Reynolds number considered 
is 10 - 100. 


Blockage ratio 


Blockage ratio defined as the flow passage blocked due to the presence of the solid 
obstacle relative to the reactor tube radius. It clearly gives us an idea how the 
substrate size, compared to the reactor dimension, affects the flow pattern and, 
in turn, the uniformity of the deposited material. 


Block position 


In the actual reactor the substrate should be placed such that the reactants react 
in front of it and the product species get deposited on the desired surface. If 
it is too close to the injector nozzle plane, this can lead to reactions behind the 
substrate surface. Too large a distance from the nozzle leads to an excessive loss 
of the product material before it is deposited on the substrate. 


Velocity ratio of the central and peripheral jets 


For the mixing of the reactant species at a favorable position in relation to the 
position of the substrate the jet velocity ratio should be properly selected for a 
certain geometric configuration. This ratio also dictates the availability of the 
reactant species to each other in the region before the deposition surface. 


Jet angle 


The angle at which fluid is injected from the central and the peripheral jet is quite 
important in the sense that it significantly influences the extent of mixing. It is 
expected that if the peripheral jet is injected towards the central jet mixing of 
the streams occurs near the inlet plane. But if the jets are placed parallel to each 
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other then the mixing and in turn the reaction zone is shifted to a downstream 
location. 


2.5 Shear stress on the solid surfaces 

Wall shear stresses in internal flows indicate a scaled magnitude of the 
velocity gradient at the solid surface. It is likely that the uniform deposition of 
the product species will be strongly correlated with a uniformity in shear stress 
at the substrate surface. This is because changes in the velocity gradient will 
also affect the rate of deposition. Moreover at high levels of the shear stress, 
the species particles do not get the required residence time to deposit uniformly 
leading to a nonuniform layer of the deposited material. Thus shear stresses 
at both the tube walls and the substrate surface are important for accurately 
predicting the flow patterns and the deposition characteristics. 
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Chapter 3 


Numerical Solution 


In this chapter numerical techniques utilized to solve the complete set of 
equations for flow and mass transfer are discussed in detail. The first section 
is devoted to orthogonal grid generation that is important for computation in 
complex geometries. The finite volume method for solving the Navier-Stokes 
equations is described subsequently in detail. 


3.1 Orthogonal Grid Generation 


The success of numerical simulation of a system of partial differential equa- 
tions is strongly dependent on the quality of the grid. For numerical simulation of 
flow in complex geometries, body conforming curvilinear grid lines are necessary. 
This is essential to accurately satisfy the boundary conditions. The exactness 
with which the boundary conditions are satisfied significantly influences the nu- 
merical solution. There are well-established techniques available in the literature 
Chikhliwala and Yortsos (1985), Duraiswami and prosperetti (1992), Nair and 
Sengupta (1998) and Beale(1999) to generate body conforming coordinate sys- 
tem. Examples are conformal mapping, transfinite interpolations and methods 
based on elliptic PDEs. In the present work we have used an elliptic partial 
differential equation for grid generation. 
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3.1.1 Elliptic Grid Generation 


The generation of boundary-conforming coordinate system is accomplished 
by determination of the curvilinear coordinates in the interior of a physical do- 
main from specified values (and/or slope of the coordinate lines at the boundary) 
on the boundary of the domain; one coordinate remains constant along each 
curvilinear line while the other varies monotonically on that line. 

The problem can be seen (in 2-D) as finding the locations of' “field points” in 
the solution domain, which correspond to regularly spaced grid points in a rect- 
angular transformed (^,?7) geometry. The generation of field values of a function 
from boundary values can be done in various ways. One approach is to use inter- 
polation between the boundaries. Alternatively, the field values can be thought 
of as solution of a boundary-value problem. The essence of this technique is to 
solve a set of elliptic partial differential equation with complete boundary corre- 
spondence. Essentially, when coordinates on all the boundaries are completely 
specified such as 0 < .^ < 1, 77 = 0, 1 and 0 < rj < 1, = 0, 1, a well-posed 

Dirichlet problem is obtained. Some of the important characteristics of the ellip- 
tic system with Dirichlet boundary conditions are listed below 


1. The system of elliptic PDEs obeys the extremum principle Muralidhar and 
Sundararajan (1995), stated as follows : 


Let the governing equation be of the form 

= P over D 


(3.1) 


where is the Laplacian operator and P is a prescribed function of x and y. 


(a) If P > 0, the maximum values of ^ occur on the boundary dQ ; 


(b) If P < 0, the minimum values of ^ occur on the boundary dQ ; and 
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(c) If P = 0, both the maximum and minimum values of ^ occur on the boundary 

aa. 


Property (c) is desirable as it ensures the generated grids represent one-to-one 
mapping. If the linear operator chosen for grid generation violates the extremum 
principle, curvilinear lines may fold back into the physical domain. This is clearly 
not acceptable because the extrema that lie inside the domain will lead to an 
overlap of the curvilinear coordinate lines. 


2. Elliptic operator has an inherent smoothness property that prevails over the 
solution derived for the whole domain. Specifically boundary discontinuities do 
not propagate into the domain. 

3.1.2 Orthogonal grid generation 

The elliptic generation system generates boundary conforming coordinate 
lines but does not ensure orthogonality of the grid lines throughout the solution 
domain. At the boundaries the grid line slope can be however specified as or- 
thogonal. Orthogonality is an important property of the curvilinear coordinate 
system required from the computational point of view. By maintaining orthogo- 
nality of the grid lines the following advantages can be achieved Ryskin and Leal 
(1983), Thompson et al. (1985), Eca(1996): 


1. Grid distortion can be minimized which in turn reduces the possibility of 
incurring greater truncation error in the finite difference expressions. 


2. The same order of accuracy for all the terms in the governing partial differential 
equations can be maintained to the level set by the discretization techniques 
associated with the finite volume formulation. 


3. With the orthogonality constraint the grid generation equation has fewer terms 
and thus reduces the amount of computations required. 
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The derivation of the orthogonal grid generation equations are outlined in 
appendix C. The two dimensional form of the partial differential equation used 
in the present work is written as : 


d ,j.dx d 1 dx 

A ( i ^)=0 

ae’ de,’^ dn^fdr,’ 


(3.2) 

(3.3) 


where / is the distortion function given by 


(3.4) 


The solution x{^, rj) and rj) of these PDEs provide the mapping between the 
physical {x, y) domain and the transformed (^, rj) rectangular domain. The two 
PDEs are supplemented by boundary conditions that assign 77 ) values to the 
boundary points (xb^yb)- 

3.1.3 Numerical solution 

Equations ( 3 . 2 ) and (3.3) are strongly nonlinear and can be solved by an 
iterative technique in the rectangular transformed domain. In the present study 
the generation equations are solved by the finite difference technique with all 
the derivatives being discretized by second order difference formulas. The most 
important issue in the solution of these equations is the determination of the 
distortion function /(^, 77 ) in the solution domain. There are three possibilities 
for calculating the distortion function and they are listed below Eca(1996): 


1. Evaluate / from its defini|;ion expression (3.4) at all the points in the domain. 


2 . Determine / from its definition expression over the boundaries where the 
physical coordinates are completely given. Next evaluate / in the solution domain 
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by a suitable algebraic interpolation scheme or by solving a Laplace equation for 
/ given by 

vV = o 


3. Choose / as a suitable algebraic function. This approach calls for experience 
and skill as it crucially depends on the geometry and the mapping layout. 

In the present study we have adopted the second alternative as this ap- 
proach involves lesser amount of calculations. As equations (3.2) and (3.3) are 
nonlinear, they are linearized by evaluating the distortion function / by values 
of coordinates from the previous iteration. The orthogonality constraint on the 
curvilinear coordinates are imposed by making the off-diagonal elements of the 
metric tensors to vanish which yields the form of orthogonal boundary condition 
as 

+ ViVv = 0 (3.5) 


Clearly this boundary condition is of nonlinear type. It is satisfied iteratively 
inside the global iteration cycle that solves the generation equation. It is to be 
noted that for a vertical boundary, x coordinates are kept fixed and y coordinates 
are calculated by the relation (3.5). For a horizontal boundary, the procedure is 
reversed. 


The solution algorithm for orthogonal grid generation can now be summarized as 
follows: 


1. Develop an initial approximation for the grid. In the present work, algebraic 
interpolation has been used. It has been observed that a reasonably good initial 
guess leads to a fewer number of global iteration cycles. 


2. Determine the distortion function at the boundaries and evaluate the function 
for the whole domain by using transfinite interpolation. 


32 



3. Solve equations (3.2) and (3.3) with a fixed by an iterative technique. 

Satisfy all the boundary conditions. 

4. Return to step 2 if convergence criterion is not satisfied. 

It has been seen that if at all the boundaries, the coordinates are specified, 
the solution converges rapidly, while the orthogonality criterion at the boundaries 
impedes the rate of convergence. Secondly, it is a well-established fact that near 
a convex boundary grid lines get attracted towards each other and near a concave 
boundary they get repelled by each other. Due to attraction-repulsion between 
the grid lines, the non-uniformity of the grid lines increases significantly. This is 
detrimental to the efficient numerical computations of the governing equations. In 
the present study we solve the set of governing equations in the physical domain 
itself, and the magnitude of the metric tensors of the generated grid do not affect 
the flow field solution. Thus the grid generation equations have not been solved 
to machine accuracy. Instead, coarser convergence limits at which the grid had an 
acceptable distribution were employed, so as to get relatively undistorted grids. 


3.2 Solution of the Mass and Momentum equa- 
tions 


In the present study the momentum equations are solved by finite volume 
method in the physical domain itself. This is justified because in the transformed 
domain the governing equations become quite elaborate and a heavy computa- 
tional effort is required. We have worked with cartesian components of veloci- 
ties and the variables have been arranged over a non-staggered mesh. Solution 
in non-staggered grids can lead to pressure-velocity decoupling, in turn result- 
ing in a pressure field completely independent of the velocity Patankar (1980). 
To avoid the velocity-pressure decoupling associated with calculations in non- 
staggered grids, momentum interpolation proposed by earlier authors Rhie and 
Chow (1983), Majumdar (1988), Miller and Schmidt (1988) and Choi et al. (1993) 
has been used. The discretization of the various terms and the solution strategy 
is discussed in the following sections. 
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3.2.1 Governing equations in integral form 


All the transport equations we encounter can be cast as the generic advection- 
diffusion equation with a passive scalar being transported in the fluid domain. 
This scalar is an intensive property, for example momentum per unit mass, en- 
thalpy and concentration. The Navier Stokes equation in the laminar, incom- 
pressible flow regime for an arbitrary control volume in the physical domain can 
be written in the integrated form of a generic advection-diffusion equation as 
follows : 

— I pdV -1- / pu.dS = 0 (3.6) 

di Jcv Jcs 


dt 


'cv 


p(j)dV + / {pu(j) — r<^V(?i).dS 


I cs 


S^dV 


'cv 


(3.7) 


where p represents fluid density, (j) stands for any intensive property being trans- 
ported and is the volumetric source term. Further, u is the velocity vector 
and r,ji is the diffusion coefficient of the scalar (j>. For the incompressible flow the 
above equations take the form 

pu.dS = 0 (3.8) 


dt 



- ^V(l>).dS = - 
P P 


S^dV 


(3.9) 


In present calculations for the three momentum equations u,v,w take the place 
of (p. 


3.2.2 Description of the finite volume method 

Equations (3.8) and (3.9) are evaluated over small volumes called finite vol- 
umes to generate the discretized form of the governing equations. The solution 
domain is divided into a number of contiguous hexahedral flnite volumes. The 
volumes are defined by the coordinates of their eight vertices. In complex geome- 
tries the finite volumes are formed by joining the eight vertices by straight lines. 
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Figure 3.1: A typical finite volume employed in the computation. 

A typical finite volume is shown in Figure (3.1). The coordinates of the vertices 
of the control volumes are generated by grid generation technique of section 3.1. 
The six neighbouring control volumes of P are denoted by E, W, N, S, T, B for 
the east, west, north, south, top and bottom control volumes respectively. The 
cell face centres are denoted by e, w, n, s, t, b in the same order. Edge centres are 
denoted by te, tn, tw, ts, be, bn, bw, bs, ne, nw, sw, se. In Figure (3.2a) and Fig- 
ure (3.2b) edge centres and the face centres used in the computation are depicted 
in two different planes. For computations we need values of the surface area of 
all the cell faces and the individual volumes of the cells. They are calculated as 
follows: 


Surface area 
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X 



z 


Figure 3.2: A typical finite volume cell (shaded) and its neighbours in the x — y 
and y — z planes. 

The finite volume vertices are numbered from 1 to 8 as in the Figure (3.1). The 
outward normal of a face is calculated by taking the cross product of the diagonal 
vectors of that face. Diagonal vectors are given by the position vectors of the two 
opposite vertices in a cell face. If Tj and Tj are the position vectors of two opposite 
vertices of a cell face then the diagonal vector is given by 

Tij = V^-Tj (3.10) 


Thus we can calculate the outward surface normal for all the faces of a control 
volume as follows: 


West face = 5(r42 x rsi) 
South face = |(r 46 x rga) 


East face Sg = Krsr x rse) 
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North face Sn = |(ri7 x ^ 23 ) 


Top face Sj = |(r 62 x rya) 
Bottom face Sb = |(r 5 i x r48) 


Volume: 


The volume of a cell is calculated from the coordinates of the cell vertices assuming 
that these vertices are joined by straight lines to form the six faces. It is given 
by 


3r7i.(S. + S, + Si,) 


The grids used for the calculation of the governing equations are shown in 
Figures (3.3)-(3.5). Three levels of grids are shown. For the reactor geometry 
with the curved substrate described in Chapter 1, partially converged as well 
as fully converged results are shown. From Figures (3.4) and (3.5) the issues of 
orthogonality and attraction-repulsion of the grid lines can be demonstrated as 
follows. When the elliptic equations for generating orthogonal grids (equations 
3.2 and 3.3) in section 3.1 are solved to machine accuracy, grid lines become or- 
thogonal to one another but the attraction near a convex boundary and repulsion 
near a concave boundary lead to severe non-uniformity of the grids as shown in 
Figure 3.5. Hence in the present study these equations are solved upto partial 
convergence and the results are shown in Figure 3.4. This is acceptable in the 
sense that non-uniformity is not encountered and grid lines are near-orthogonal. 
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(a) 



x/D 


(b) 



x/D 



x/D 


Figure 3.3; Grids used in the calculations for the reactor with plane substrate; 
(a) 73 X 38, (b) 83 x 48, (c) 103 x 68 
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(a) 



x/D 


(b) 



x/D 


(c) 



x/D 


Figure 3.5: Fully converged grids for the reactor with curved substrate; (a) 73x38, 
(b) 83 X 48, (c) 103 x 68 
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3.2.3 Discretization Procedure 


Here we present the discretization of various integrals that appear in the conti- 
nuity and momentum equations (3.8, 3.9). The treatment adopted here follows 
closely the report of Satya Prakash et al. (1996). 

Continuity Equation 

The integral form of the continuity equation is given by 

pu.dS ft; ^p(u.S)j = J]p(uj.Sj) = (3.11) 

J 3 j 

where Fj is the mass flux through the face j, defined by 

Fj = pUj.Sj (3.12) 



Thus the discretized continuity equation states that total mass flux through all 
the cell faces vanishes for a finite volume. In turn the net fluid mass in a control 
volume remains unchanged, a statement of the principle of conservation of mass. 
The discretized version of the continuity equation is given by 

y ^ Fj = Ejy + Fs Fg + Fn + Ft + F), = 0 (3.13) 

3 

Momentum Equations 

The momentum equation when written in the generic advection-diffusion form, 
consists of four terms; unsteady, convective, diffusive and the source integral 
related to the pressure gradient. Individually, the discretization of the four terms 
are performed as follows: 


Unsteady term 
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The unsteady term in the governing equation is discretized with the help of the 
basic assumption that all the dependent variables for which we seek solution are 
defined at the cell centroid. Thus it follows 


dt 


p(j)dV 




At 


pVp 


- (j)'^ 
At 


(3.14) 


where Vp is the volume of a typical control volume and At the time increment 
employed for the time marching. The time discretization of equation (3.14) has 
an error of first order in time. 


Convective term 


The convective integral in the governing equation is approximated as following 

pu^.dS ^(pu.S)j0j = Fj(pj (3.15) 

J 3 

where (p^ is the value of the dependent variable at the cell face and Fj is the 
mass flux through the face j. The value of (p at the cell face centre is assumed to 
prevail over the face in question. Thus the discretized equation takes the form 

pUi^.dS ~ Fyj(p^ + Fg^s + Fg^'e T Fn(pn + F^Pt + Fi(p\) (3.16) 



It is evident from the above expression that for the discretization of the convec- 
tive flux we need the meiss flux at the cell faces. This step requires velocities at 
the cell faces. It can be shown that if the mass flux at any cell face is calculated 
using velocities obtained by linear interpolation from the two neighbouring cells 
adjacent to the face, an arbitrary velocity field with a checker-board distribution 
may result, a completely nonphysical situation Patankar (1980). The reason for 
the error is the following: When velocities at a cell face are calculated by linear 
interpolation between two adjacent cell centre values, the continuity equation be- 
comes the difference between two alternate cell centre velocities and not between 
adjacent ones. This leads to decoupling of the velocity at the cell for which the 
continuity equation is discretized and any zig-zag velocity distribution satisfies 
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the continuity equation. To avoid such velocity wiggles, the momentum inter- 
polation idea has been proposed in the literature. This technique uses pressure 
weighted interpolation of velocities at the cell faces and uses these velocities to 
calculate the mass flux at the cell face. 

If J and J-l- 1 are two neighbouring cells adjacent to a face then the mass 
flux at the face j is calculated by velocities obtained by straightforward linear 
interpolation as 

Fj = pU,.S, ( 3 . 17 ) 


where Uj is the velocity linearly interpolated from the two adjacent cells to the 
face J, written as 


TT = 

' + Fj+i 


(3.18) 


The momentum interpolation adds a pressure gradient term with the linearly 
interpolated part of equation (3.17), namely 

Fj = pUj.Sj — AtVp.Sj (3.19) 


This technique serves two purposes as a whole. First it excludes the possibility of 
any velocity wiggles because velocities at a cell face depends on the pressure at 
the adjoint cell centres which in turn depend on the flux at the cell faces for which 
continuity equation is discretized. So any abrupt changes in the velocity field can 
be sensed by the continuity equation and an arbitrary velocity field cannot be 
permitted as the solution. Secondly the pressure gradient term which takes val- 
ues from the two adjacent cells (J and J -1- 1) behaves like a natural force driving 
the mass flow at the cell face j, a key point for the success of the staggered grid 
arrangement. 

As the velocities are defined at the cell centroid we have to use an inter- 
polation scheme to obtain velocities at the cell face centre. In the present work 
this has been accomplished by taking contribution of both second order central 
difference and first order upwind difference scheme. The central difference scheme 
can lead to numerical instability and the upwind method introduces numerical 
diffusion, when they are individually used. Thus using the right combination of 
these two schemes is expected to improve the robustness of the algorithm. For 
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strongly nonlinear and convection dominated flows the upwind method maintains 
the transportive property of the flow, and the central difference scheme being of 
second order reduces error in the solution. Here we demonstrate the complete 
form of the two schemes and their joint application in a finite volume code. The 
central difference method applied to calculate the face centre value of the depen- 
dent variable leads to 


(j)j — 


Vr 


Vj + 16+1 


<l>J+l + 


V 


.7+1 


6 + V>+i 


<t>J 


( 3 . 20 ) 


where j and J represent the face centre and cell centre values. Here J and J -t - 1 
are two neighbouring cells adjacent to the face j. It is to be noted that for uniform 
grids the above formula is second order accurate but for non-uniform grids this 
volume weighted formula gives only first order accuracy. The upwind scheme 
rightly manifests the transportive property of the fluid flow equations in the local 
flow direction. It introduces numerical diffusion in the solution and hence leads 
to smearing of a sharp front. This scheme is written as follows: 

F,4>j = \F„0\<t>J + \-Fj,0\cl>j+i (3.21) 


With this scheme, the dependent variable at the face centre takes the value of 
the upstream cell if the mass flux is positive and takes the value of downstream 
cell if the mass flux is negative. The two schemes on combination take the form 

+ I - F,.0|j.jti)+7[q( ,. 6/ 

Vj + Vj+i \'j -t- Vj+i 

(|f,.O|0j+l-F„O|.^j+i)) (3.22) 


where the factor 7 signifies the contribution of each scheme. It varies between 
0 and 1 ; the limit of 0 implies a purely upwind and 1 implies the purely central 
difference scheme. There is no universal rule to choose this factor in a particular 
calculation. It depends on the nature of the problem, numerical technique being 
used and the range of parameters taken into account. In the present study we 
have used 7 =0.5. 

Using equation ( 3 . 22 ) the convective flux at all the cell faces are calculated 
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as follows 


Fu)4’w — 01</>H^ — 1 — Fm, 0|</>p) + 7 [Fu,(— — ■ " (^p + — — ^ - (f)v/) 

Vw 'T yp Vw + I'p 

“(IFi,, Q\4>w — I ~ Fu, oi</>p)] 


F0S = 


(1F„ Ol-^p - I - F, 0|(/.s) + 7[F(^7-^</)p + 

Vs + Vp 

-(|F,0|<^p-1-F,0|<^5)] 


v> 

Vs + Fp 


(l^s) 


Fe<i>e 


(If., Ol^p - 1 - f., 0|« + 

-(lf„0|,^P-|-f.,0|fc)] 


Fn 4 ^n 


(|f„, o\<Pn - 1 - f„, o|« + -y[Up^<t‘p + ,, M 
-{\Fn,0\<l>N-\-FnA^p)] 


Ft<pt = (IF, 0|(?i>p - I - F, Ol(/>r) + 

-(lF,0 l<^p-l-F,0 |,/.r)] 


Vn Vp 

= (|Fb,0|.^B - I - F,0|^p) ^ y^ 4>B) 

-{\Fb,mB-\-FbM4>p)] 

Diffusion term 

The diffusion flux of the variable (f> in the governing equation can be ap- 
proximated as 

[ r^V(i>.ds « = Y. Ff (3.23) 

j d 
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where the superscript d denotes diffusive flux. Discretization of this term in com- 
plex geometries with curvilinear grids is carried out as follows. The diffusion 
mechanism is considered in the normal and tangential directions with respect 
to the curvilinear coordinate lines. These are termed normal and cross diffu- 
sion. For any face the outward drawn normal vector can be written as the linear 
combination of three linearly independent unit vectors which are themselves not 
necessarily orthogonal. Thus 

Sj = ftiHi + 01.2^2 + QjsHs ( 3 - 24 ) 


where ni, n2, 113 are three linearly independent unit vectors given in terms of their 
cartesian components as 

m = (mi ni2 nis) 
n2 = (mi n22 ms) 
ns = (mi m2 ^33) 


Thus the vector product of equation ( 3 . 23 ) can be written as 

= V^.(Q;ini -h Q; 2 n 2 + 0:3^3) 

= Q!i V(^.n.i -|- q; 2 V<5ii.n2 + 0:3 V(^.n3 ( 3 . 25 ) 


Now if A<^i, A(/i2, A</»3 are the differences in <j) along the three line segments 
Axi, Ax2, Ax 3 then 

A01 = V(j).Axi, A^2 = V(^.Ax2, A03 = V(?!>.Ax3 


If Axi,Ax2,Ax 3 are in the directions of ni,n2,n3 and Axi, Aa:2, Aa;3 are the 
magnitude of the line segments respectively then 


A(l)i = V(^.niAxi, A4>2 = V^.n2Aa;2, A(^3 = V(j).nzAx3 


Thus it follows that 


V(^.ni 


A01 

Axi’ 


V(f>.n2 


A^2 

AX2' 


V^.ns = 


A 03 

Axz 


( 3 . 26 ) 
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The diffusion flux at the face j can be written using equations (3.25) and (3 26) 
as 

T = = r,(a. II + (3,27) 

Equation (3.24) can be written in matrix form as 


nil 

n2\ 

nz\ 


5i, 

ni2 

n22 

nz2 

a'2 

= 52, 

. niz 

'/^23 

nzz _ 

1 «3 J 

1 1 53 J 


where Sij,S 2 ],Szj are the three cartesian components of the outward normal 
surface vector Sj. The coefficient matrix in the above equation can be inverted 
to get the values of Q;i,Q; 2 )Ct 3 - We have used Cramer’s rule to invert the matrix. 
Hence 


ai 


£i 

D ’ 


a'2 


D ’ 


as = 


D 


(3.28) 


where D is the determinant of the coefficient matrix of the matrix equation. The 
matrices Di, D 2 and D 3 are formed by replacing the respective columns of the 
coefficient matrix by the right hand side vector of the matrix equation. 

For a system of linear equations 


Ax = 6 


with n unknowns, if the rank of the coefficient matrix is equal to n written by 

rank A = n 


then the coefficient matrix A is non-singular and the linear system has a unique 
solution Kreyszig (1993). The row vectors of the coefficient matrix of equation 
(3.24) are linearly independent. Thus the maximum number of linearly indepen- 
dent row vectors are three, resulting in the rank of the matrix to be three. Hence 
this matrix is non-singular and gives a unique solution for the a’s. 
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Figure 3.6: Lengths required for the calculation of diffusion fluxes and related 
nomenclature. 

In the discretization of the diffusive flux, the normal and the cross diffusion com- 
ponents have to be accounted for. It can be shown that for orthogonal cartesian 
grids the cross diffusion term vanishes but for non-uniform curvilinear grids, cross- 
diffusion in the three directions are possible. In Figure (3.6) the diflfusion fluxes 
and the related terminology are clearly shown. 


For the calculation of the cross diffusion fluxes, values of the dependent variables 
at the edge centres are required. These values have been calculated by the volume 
weighted interpolation technique using cell centre values of four neighbouring cells 
surrounding an edge. Values at the four edge centres on the west face are written 
below: 


VTw<i>p + yw<t>T + yp4>TW + VT<t>w , _ Vsw(l>p + Vw<Ps + Vp^^sw + Vs4>W 

" Vtw-^Vw + Vp + Vt ’ ” Vsw + vw + Vp -f Vs 
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, _ Vbw4>P + yB’t'W + yp4>BW + Vw(t> 

+ 1/,. + Vu. 


B 




_ VnW^Pp + Vn4'W + ^P<pNW + Vw 


(t>M 


Calculations for the edge centre values at other faces are similar and straightfor- 
ward. 


Pressure Gradient term 


In the momentum equation, the evaluation of the pressure gradient term is 
identical to the diffusion term. Its discretized form is written below for the ith 
momentum equation after integrating over a typical finite volume. 


Vp.n,dV = Vp.n^dV + pV .n^dV 

J V J V 

= f V.{pni)dV 
J V 

= f pnj.dS 

J $ 

j 


Thus the discretized form of the pressure flux is given by 


jvp.nidV = Y^PJS^J 


(3.29) 


This term is similar to the diffusion flux of equation (3.23), and can be approxi- 
mated in an identical manner. 


3.2.4 Time integration Scheme 

The momentum equation is mathematically of mixed character. It can be 
visualized as elliptic in space and parabolic in time. The momentum equation 
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has been time integrated semi-explicitly with each term treated explicitly except 
for the pressure term which is implicitly treated. The time integration scheme 
mainly consists of two steps. In the first step which is the predictor, velocities 
are computed based on pressure field at the current time level. In the second, 
the corrector step, velocities and pressure are corrected iteratively to satisfy the 
divergence-free condition for the velocity field. The pressure-velocity coupling 
in the corrector step reduces to a Poisson equation for the pressure correction. 
This equation has been solved iteratively until the continuity equation is satisfied 
to the machine accuracy. It should be noted that an explicit formulation does 
not involve any matrix inversion for the velocity calculations compared to semi 
implicit type of formulations such as SIMPLE and SIMPLER Patankar (1980). 
Thus explicit calculations are not computationally intensive. However the explicit 
formulation suffers from the time step restriction due to its conditionally stable 
behaviour. On the other hand, very small time step increment can resolve the 
minute details of the transients. Weak coupling between discrete points in the 
domain is however inferior compared to an implicit formulation. Although the 
present study has not been directed towards obtaining the transient details, the 
points mentioned above justify the choice of the explicit algorithm. The predictor 
and the corrector steps with the equations to be solved are written below. 


Predictor 


In the predictor step, the velocities (marked by *) are predicted using values from 
current (nth) time step: 

J 3 

Equation (3.30) leads to the predicted velocities as 

u-p = u"p - (3.31) 


Corrector 



(3.32) 


By advancing semi-explicitly in time the momentum equation becomes 

+ Epr'sj 

J 3 

By subtracting equation (3.31) from equation (3.32) we get 

= (3 33) 

where u'p and p' are the velocity and pressure corrections to the predicted veloc- 
ities and guessed pressure. They are defined as 

u'p = - u), (3.34) 

p' = (3 35) 


If equation (3.33) is cast into the integral form then a relation between velocity 
and pressure corrections can be obtained as follows: 



with the integration being performed over an arbitrary finite volume. The inte- 
grands of each side can be pulled out of the integral to get 

u' = -— Vp' (3.36) 

P 

In the corrector step we enforce continuity implicitly, written symbolically as 

V = 0 ^ (3.37) 
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Writing this term as sum of predicted fluxes based on predicted velocities of the 
first step and the amount of flux to be corrected to get a divergence free flow 
field, we get 

'£p; + T,F', = C' (3.38) 

J 3 

This gives 

3 3 

Thus equations (3.36) and (3.37) constitute the corrector step. These two equa- 
tions are to be simultaneously solved. It has been accomplished by an iterative 
technique. The two equations are reduced to a Poisson equation in pressure 
correction, solved to get a divergence free velocity field. 


Pressure Correction equation 


It can be shown the two equations to be solved in the correcter step (equa- 
tions 3.36 and 3.37) lead to a Poisson equation in pressure correction. The two 
terms of equation (3.38) when written in integral form lead to 




pu*.dS 

fv.a’dV 


= -f pu'.dS 

j Jcs 


= -p 

[ V.n'dV 


lev 

= At 

[ ^.{Vp')dV 


lev 

= At 

[ V^dV 


lev 
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Thus the integral form of the equation (3.39) is 


At f V^dV = p 

J V 


V.n*dV 


As the integration has been carried out over an arbitrary finite volume, the inte- 
grand of each side can be pulled out of the integral to get 

VV = ^V.u* (3.40) 


Hence solving equations (3.36) and (3.37) is equivalent to solving equation (3.40), 
the Poisson equation in pressure correction. The two equations to be solved 
simultaneously in the corrector step can also be reduced to a form which is more 
amenable to the Gauss-Seidel iterative method. This is called the residual form 
of the pressure correction equation and is developed as follows. 

The mass flux correction can be written in terms of pressure correction by 
using equation (3.36) as 


F'j = pu'j.Sj 

= —AtVp'.Sj 


It follows that 

^F', = -5]AtVp'.S, (3.41) 

3 3 

The pressure correction term in the above equation can be discretized in the same 
way as the diflfusion term: 

Ai Y, VP'.S, = A(l5J^(Vp'.S,)„j + 53(Vp'.S,).i] (3.42) 

3 3 3 

where the subscript nd and cd corresponds to normal and cross diffusion terms, 
explained in the context of the diffusion flux. It is to be noted that the cross 
diffusion term contains values of pressure correction from both the current as 
well as the previous iteration levels. The normal diffusion term is discretized in 
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the same way as the normal diffusion flux. Using equation (3.27) we write 


At = At[ 

j 

Qf2 


0^1 

Aa:i 


(p'p-pV) + ^ Is (p's-p'p) + t^ le (p'e-pV) 


ai 


Ax 


Axi 




Q3 

Aa-3 


Thus equation (3.42) can be written as 

Ai 52 Vp'.S, = Ai ^ ^ p'„, + appV + Air, 


n5 


cd 


(3.43) 


where nb denotes all the neighbours of the cell with centroid at P. The coefficient 
ap is given by 

ap = At(-^ 


1 1 

ai 

. ^2 1 

0^3 1 

«3 , 

Ihl A 5 

Ax2 

Axi 

‘ A ^ 

Aa:2 

Axz ‘ 

AX3 ' 


Residual form of the pressure correction equation 


In iterative techniques often it is useful and advantageous to solve the difference 
equation in residual form Tannehill et al. (1997). This is demonstrated below. 

Any difference equation for a variable can be written as the sum of 
the difference equation of a provisional quantity and an incremental quantity 

L 4>i,, = L + L 64>^J = 0 (3.45) 


The provisional solution 4>i^j represents a value that might occur at some point 
in an iterative process. The correction is the quantity to be added to the 
provisional to get the full solution. The deflnition of the correction term may bear 
a slightly different meaning, depending upon the algorithm and the nature of the 
problem. These two quantities give the exact numerical solution when combined. 
The residual R is defined as the number obtained when the difference equation, 
written in a form producing zero on the right hand side, is evaluated on the basis 
of the provisional solution. This is given by 

R = L\^ (3-46) 
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In the corrector step the mass flux at any face j at the (n+l)th time step summed 
over all the faces of a finite volume is written in an operator form as 

= L = 0 (3.47) 

3 

Using equation (3.45) we get 

L = LFj + L5F^ (3.48) 

The term, Fj denotes the provisional flux, namely the quantity at a certain itera- 
tion level. The symbol 5Fj is the correction term to be added with the provisional 
solution to reach the exact numerical solution. The provisional quantity can be 
written in terms of the predicted value at the start of the iteration (which does not 
change during the iterations) and the quantity that has modified the predicted 
flux to the provisional one during the iterations already occurred. Hence 

= LF^ + L6Fj 

= lf; + l (f;)^ + l sf, ( 3 . 49 ) 

with q representing the current iteration level. According to the definition of the 
residual given by equation (3.46) above, equation (3.49) gives the value of the 
residual as 

r = -lf;-l (F,)^ (3.50) 

Thus by equations (3.47), (3.49) and (3.50) 

L8F, = -lf;~l{f;)^ 

= R (3.51) 

It is to be noted from equation (3.48) that as the provisional solution reaches the 
exact numerical solution F, the residual vanishes. 

The incremental flux 6Fj is the difference between the values at the next 
iteration level and the current iteration level. Using equation (3.43) this term 
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can be deduced in terms of pressure corrections at the cell of interest and the 
neighbouring cells. If pressure corrections at all the neighbouring points are 
neglected then the incremental flux is given in terms of the pressure correction 
values at point P at the current and next iteration levels as follows: 

L 6F, = a,[(Pp),+, - {Pp),l (3.32) 


It follows from equations (3.51) and (3.52) that 

R (3.53) 


Equation (3.53) written with the application of a relaxation factor takes the form 


(Pp)9+i 


{P'p)g + 7 


A 

Gp 


(3.54) 


Equation (3.53) is the residual form of the equations (3.36) and (3.37) to be solved 
in the corrector step. It can be used to update the pressure correction during the 
iterations for getting a divergence-free velocity field within a time step 

During the derivation of the residual form, pressure correction at all the 
neighbouring cells have been neglected. However this is just an iterative strate- 
gies ti obtain corrections and the field solution will enforce continuity of the entire 
discretized equation (3.39) without neglecting any terms as the full equation is 
used to compute R. This is in principle similar to the exclusion of all the neigh- 
bouring velocity corrections in the SIMPLE algorithm Patankar (1980). This is 
justified because with the inclusion of pressure corrections at all the neighbour- 
ing cells pressure correction updating expression becomes unmanageable. As a 
neighbouring cell will bring the influence of its neighbours and its neighbours will 
bring the influence of their neighbours and so on. The relaxation factor mentioned 
above is applied to accelerate the rate of convergence of the pressure correction 
equation. This factor may bear a significant meaning depending upon the nature 
of the system of linear equations arising from the pressure Poisson equation. 


The solution algorithm for the momentum equations is now summarized. 
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Initial conditions for velocities are given and pressure is guessed for all the points 
in the domain to start the calculation. 

Obtain velocities and the mass flux in the predictor step using the initial condition 
or values from the previous time step as; 
start: 1 For each cell (P) 

3 3 


close: 1 

start: 2 For each face {j) of a cell (P) 


Pj = pUj.Sj — AiVp.Sj 


(6j is the linearly interpolated velocity at the cell face given by equation 3.18) 
close: 2 

Boundary conditions for the velocities and pressure are applied. 

Initial guessed values for pressure correction are utilizes to start the iteration in 

the corrector step. 

start: 3 Start iterations {k): 

Calculate corrective mass flux F'^ for all the cell faces, the residual R at each cell 
and pressure correction values are updated as 
start: 4 For a face (j) of a cell (P) 

p; = -Atvp'.Sj 

close: 4 

start: 5 For a cell (P) 

3 3 

close: 5 

start; 6 For a cell (P) 

, , R 

p •<- p + 7 — 
ap 


57 



close: 6 

check for convergence as 
If 


(where e is the prescribed convergence criterion, say 10~® %) 

STOP 

Else: Continue 
close: 3 {k) 

Update velocity and pressure 
start: 7 For a cell (P) 

u u + u' 

p -f- p + p' 

close: 7 

Repeat the same calculations for each time step till steady state is reached 

3.2.5 Calculation of the time step 

It can be shown by the von Neumann stability analysis that for a three di- 
mensional advection-diffusion equation with an explicit FTCS (Forward in Time 
and Central in Space) scheme, the numerical error does not get amplified in suc- 
cessive time steps if the following two conditions are satisfied : 

A. CFL criterion 



It states that fluid flow information should not propagate by a distance greater 
than one cell in a given time step. This is the convective limitation which the 
numerical algorithm should obey to be stable. It is written as 
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Here u, v and w are the three components of velocity along the three cartesian 
directions x, y, z respectively. If we write this criterion in terms of nondimensional 
variables it remains unchanged with the nondimensional quantities replacing the 
dimensional quantities. 


B. Grid Fourier number criterion 


Fourier numbers based on individual grid size should not exceed a certain value 
for the numerical scheme to be stable. This criterion is wiitten as: 


A / 1 1 1 N 1 


with the variables as stated above. If this criterion is written in terms of the 
nondimensional variables it takes the form 


At 1 1 

Re^A2:2 + ^y2 ^ A22'' - 3 


(3 56 ) 


More generally, the kinematic viscosity u is replaced by the diffusivity (thermal 
or mass) and Re by Pe, the Peclet number. The above two criteria leads to a 
single expression as 


At < 


2 




( 3 . 57 ) 


For calculations over non-uniform grids where cell lengths in the three directions 
vary through out the domain, dimensions depicted in the above criteria have to be 
calculated by including all the cells. In the present study we have calculated the 
minimum of the six diffusion lengths for the six faces along a cartesian direction 
and then obtained the minimum for that cell. The minimum of the diffusion 
lengths for a cell is then calculated by comparing the lengths of all the faces. 
Finally the minimum of the diffusion lengths for the whole domain is calculated. 
Similar calculations are carried out for the other two cartesian directions. These 
three lengths have been taken as the dimensions that dictate the propagation 
of information in the flow domain. The following expression for calculating the 
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time step has been used after obtaining the appropriate lengths that govern wave 
propagation. 


At = w 


1 


( Ji I V 1 I 3_/_l I 1 

' Ai Ay ~ Az^ 1 Ai^ ' Ay^ 



(3.58) 


Here w (< 1), is a factor that makes the time step more conservative than what is 
demanded by the stability limit. It has been observed that the above expression 
works well for all Reynolds numbers for stability factors of the order of w = 0.9. 


3.2.6 Initial and Boundary conditions 

As our aim is to obtain a steady state, we can start with an arbitrary initial 
condition. However to obtain a meaningful transient variation, initial condition 
should be carefully selected in accordance with the physical situation. 

As all the dependent variables are defined at the cell centres and the bound- 
aries lie on the cell faces, a fictitious layer of cells adjacent to each boundary 
segment have been employed in order to apply the boundary conditions. By in- 
terpolating values between a fictitious cell and a real cell adjacent to a boundary 
the values at the corresponding boundary can be calculated. It is important to 
note that if the fictitious layer of cells are perfect images of the real cells adjacent 
to the physical boundary, the boundary conditions are exactly satisfied. Due to 
the curvilinear non-uniform grids in the present study volume-weighted interpo- 
lation has been used to set all the boundary conditions. Boundary conditions 
applicable to all the boundaries have been discussed in Chapter 2, and will not 
be repeated here. Only the boundaries where special treatment is required are 
discussed below. 


Corner cell treatment 


The two corner cells of the solid obstacle are different from the rest of the bound- 
aries in the sense that they have more than one face exposed to the fluid flow. In 
the present work we have set no mass penetration conditions at these faces. This 
condition is outlined for one cell as in Figure 3.7 
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Figure 3 . 7 ; One of the corner cells (shaded) and associated velocities. 


The normal velocity at the face b is calculated by using volume weighted 
interpolation of velocities between the cells 1 and 2 . Setting this velocity equals 
to zero leads to 


uicoso;! + uisinQ:i = — ^(uscosoii + U2SinQ'i) 

V2 


In the same way normal velocity at the face a can be calculated and when this is 
set to zero we get the following relation 


—ui sin 0:2 + cos 0:2 



{-U2 sin 0:2 + V3 cos 0:2) 


Solving the above two equations for ui and Ui we get the expressions for the 
velocities at the corner fictitious cells as 


Ui 


— ^ (us sin Oil sin 0:2 ~ '^^3 sin cki cos 0:2) — ^ (^2 cos 0:1 cos 0:2 + ^2 sin ai cos 0:2) 


COs(Q!i - 0:2) 


( 3 . 59 ) 


61 



Vi 


I/3 cos Oil sin ct2 — V3 cos cti cos 02) — ^ (u2 cos ai sin 012 + V2 sin ai sin 0:2) 


cos(q!i — 0^2) 


( 3 . 60 ) 


Periodicity boundary conditions 



Figure 3 . 8 : Velocities associated with the free-slip boundary condition in the 
azimuthal direction. 

As the computational algorithm of the present work has been generated for 
the solution of a three dimensional advection-diffusion equation, special treat- 
ment is required to make it applicable for a two dimensional geometry. In the 
present study the free-slip boundary condition has been applied in the third az- 
imuthal direction. This condition states that the normal velocity and normal 
gradient of all other variables including the tangential velocities should vanish. 
Implementation of this condition is demonstrated with the help of Figure 3.8 for 
the left face. 
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Normal and tangential velocities are calculated for both the real and ficti- 
tious cells using the discretized azimuthal angle ^ per cell as 

V” = —Vi sin ^ H- Wi cos ^ 

V 2 = —V 2 sin ^ + W 2 cos - 

-tr ^ Txr ■ ^ 

— Vj cos “ I'V 1 sin — 

j’i T/ ^ Txr ■ ^ 

Vo — 1'^ cos — -l- 14'^2 sin — 

^ 2 2 

Using these velocities the normal component of velocity and the normal gradient 
of the tangential velocities can be calculated. Setting these two quantities equal 
to zero, the boundary conditions for the velocities can be obtained as follows: 

Ui = U2Cos0 + W2sin^ (3.61) 


Wi = V 2 sin 9 -W 2 cos 9 (3.62) 


As the axial and azimuthal directions are orthogonal to each other, boundary 
conditions for the axial direction are straightforward in the sense that the axial 
component of velocities can be directly used. Boundary conditions for the right 
face in Figure 3.8 can be similarly obtained. 


3.3 Solution of species transport equations 

In Chapter 2 on mathematical formulation of the species transport equa- 
tions, both dimensional and nondimensional forms were presented. These equa- 
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tions use the mass- averaged velocity obtained by solving the Navier Stokes equa- 
tions. The equations governing the momentum transport and the species trans- 
port mechanism form an uncoupled system when their respective time scales are 
widely separated. The steady state solutions obtained from the momentum equa- 
tions can then be directly used for solving the species transport equations. For 
the sake of illustration the discretization of the transport equation for species A 
is written here in nondimensional form. 

+ V.Kpu^n - = Ra (3.63) 

at i em jja,o 

where 

Ra = - kflOAtOB) (3.64) 

is the rate of consumption of species A due to a homogeneous elementary chemical 
reaction of the form 

kf 

vaA + i^bR 7^ 

kr 

Equation (3.63) is similar to the momentum equation with the exception of the 
source term Ra- The source term is appreciable and important since it signifies 
the consumption of species A. The individual terms in equation (3.63) are dis- 
cretized in the same fashion as the momentum equations in an advection-diffusion 
framework. The discretization of the source term requires special treatment and 
is further discussed. It can be seen from the set of species transport equations 
(2.33-2.36) that they are coupled with each other by virtue of the source term. 
Moreover the source term is nonlinear and has to be linearized in order to generate 
a set of linear algebraic equations. 

3.3.1 Source term linearization 

To develop a numerical algorithm that has a traditional structure, the lin- 
earization of the source term is essential Patankar (1980). This can be accom- 
plished as follows. The source term can be approximated as being linearly depen- 
dent on the concentration, over a limited range of values. It is self-evident that it 
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is better to assign the source term to be linearly varying rather than a constant. 
This linearization formula can be written as 

R = Rc + Rp4’p (3.65) 


where (pp is the dependent variable of interest (for example, mass fraction of the 
one of species), Rc and Rp are the curve-fitting coefficients of the linearization 
step. It is required that Rp is negative. The negative-slope dependence of the 
source term on pp is illustrated by the following simple analogy. If in a system 
there is a volumetric heat source, there should be a proper heat removal mecha- 
nism. Otherwise the system will experience an indefinite increase in temperature. 
Moreover, with the positive slope variation the coefficient of the dependent vari- 
able in the algebraic equation may become negative. This is conceptually not 
acceptable when iterative schemes of matrix inversion are employed Pataukar 
(1980). A negative-slope linearization ensures that instabilities and unrealistic 
physical solutions do not arise. In the present work, the linearization is obtained 
by calculating Rc and Rp from the previous time step values. The full source 
term for species A in linearized form is written below: 


Ra = 


d PcMa^, , ^ 

^ (ferCJc^D - kfUJAUJB) 

Rc ~ Rp^A 


(3.66) 


where 


d 

Rc = ^Ajj-^A,o 


Pc Ma 
Me Mo 


k^pCOc^ D ) 


d Pc Ma, 


While discretizing the species transport equations, convective and diffusive fluxes 
are calculated explicitly and the time marching procedure is used to reach the 
steady state solution. The discretized equation is written as: 

+ x;(j?)”i + Ks - RpMT' (3.67) 

j 3 


where i denotes all the four species. The linearization coefficients are summarized 
below for all the four species transport equations: 
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species A 


p ,, ^ , 


Mr Mr 


P _ d Pc Ma, 


species B 


D _ ^ Pc Mb , d Pc Mb , 

Rc - , Rp - 


species C 

d Pc Me , d Pc Me , 

Rc = , Rp = 

species D 

_ d Pc Me, „ d Pc Md , 

It is clear from the above expressions that the slope of the linearization formula is 
negative for the reactant species since kf denoting forward reaction rate constant 
is positive for the conventions used in Chapter 2 but is zero for the product species 
for an irreversible reaction {kr = 0). 


3.3.2 Operator Splitting algorithm 

The chemical reaction we have dealt with is spontaneous and the forward 
reaction rate constant is quite large in magnitude. In the species transport equa- 
tions the source term linearly depends on this large constant. With Peclet num- 
bers associated with mass transfer being large as well 1) large value of the 
source term leads to numerical difficulties during matrix inversion. 

Though it has not been the target of the present study to resolve the reaction 
front moving towards the substrate, the mass transfer problem has been solved 
by the application of the Operator Splitting (OS) algorithm Cooke (1986), Ding 
and Liu (1989), Muralidhar et al. (1993). This technique is well-established to 
give accurate and detailed results when applied to moving front problems. While 
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most of the commonly used algorithms fail to explore the sharp variation of a 
variable in the high Peclet number region, the OS algorithm does not smear the 
fronts and estimates the transport of the front quite accurately. The basic essence 
of this algorithm is to split the governing advection-diffusion-reaction equations 
into two or three equations, one for advection, one for the diffusion and one for 
the reaction part and solve each of the the homogeneous part of the equation as 
accurately as possible. Thus the full equation which is difficult to solve in the 
presence of a large source term is split into parts which are partially solvable 
preferably by an analytical technique. Hence the splitting procedure not only 
reduces the computational effort but also ensures successful numerical simulation 
in the presence of high Peclet numbers and large source terms. In the present 
study the governing transport equations for each of the species are split into two 
parts as follows: 

^ 1 /~) 

— (pcjj) + V.[(pua;j) — — — 77 ^V(pa;i)] = 0 over a time step At (3.68) 
ot Pem Da^o 


d_ 

dt 


{pu!t) = Ri over the time step At 


(3.69) 


These two equations can also be termed as the predictor (equation 3.68) and 
corrector (equation 3.69) steps. The predictor step involves solving an advection- 
diffusion equation without any source term by marching in time, subject to all 
the boundary conditions for the multidimensional problem. The corrector step is 
solved with the initial condition derived from equation 3.68. It has a closed form 
analytical solution. The numerical treatment of the two steps are outlined below. 

The predictor step is solved by the finite volume method outlined in section 
3.2. After integrating it in an arbitrarily finite volume and approximating various 
integrals as in the momentum equations we get : 


pVp 


At 


j 3 


The time advanced value of concentration is now obtained as 


(wi)p''' 


(a;i)p - 


At 

Wp 




(3.70) 
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For the corrector step equation (3.65) is inserted into equation (3.69). By inte- 
grating equation (3.69) we get the analytical expression for the mass fraction at 
(n -1- l)th time level 

(cj,) J+i = Ai-A 2 exp(-RpAt) (3.71) 

where Ai and A 2 are given by 

= (h)p , /t2 = (3.72) 


The coefficients Rc and Rp are given by the linearization expression for each of 
the species. 

The operator-splitting algorithm for equations (3.68) and (3.69) can be sum- 
marized as follows: 


Give initial conditions for all the species in the domain to start calculation 
Obtain values of mass fraction by solving the predictor step as: 
start: 1 For a cell (P) 

Mp = Mi - 


close: 1 

Obtain values of mass fraction by solving the corrector step using values from the 
predictor step as the initial condition as: 
start: 2 For a cell (P) 


(wjp+^ = Ai - A2exp{-RpAt) 


{Ai and A 2 are given by equation (3.72) and expressions for Rc and Rp are given 
in section 3.3.1) 
close: 2 

Apply the boundary conditions. 

Repeat the calculations till steady state is reached. 

The solution from the OS algorithm was tested directly against that from 
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the finite volume code applied to the complete system of mass transfer equations. 
It was found that: 


1. The transient details from the two approaches follow each other very closely 
and the steady state distribution of species are practically identical. 


2. Computational time required is lower for the reaction-split approach by a 
factor of 1.8. 


3. With the forward reaction rate parameter kf becoming large it is difficult 
to handle all the terms of the species transport together, specifically, the finite 
volume code diverged for /c/ > 10^. The application of the operator splitting 
algorithm allowed larger source terms upto kj = 10®. In view of these advantages 
the final results of the present work for mass transfer have been obtained by the 
OS algorithm. 

3.3.3 Initial and Boundary conditions 

As the present study was targeted towards getting only the steady state 
solution, constant values have been assigned to all the mass fractions as the 
initial condition. It has been seen that when converged solution of a particular 
case is used as the initial condition for an another case whose parameters does 
not vary significantly convergence is quite faster. All the boundary conditions 
to be applied for complete solution of the species transport equations have been 
discussed in Chapter 2. Only the boundary condition at the deposition surface 
is discussed here. 

The mathematical models which consider the surface reactions, derive the 
boundary condition at the deposition surface by equating net diffusion flux at 
the surface to the rate of reaction at the surface Mahajan (1996). But for the 
present model which considers the deposition surface to be inert the specification 
of boundary condition at the deposition surface is not straight-forward. The above 
mentioned condition reduces to a no-diffusion condition for a passive surface. The 
no-diffusion condition does not pose any problem for all the species other than 
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the species to be deposited. However if this condition is also specified for the 
depositing species there is no obvious and consistent method for the calculation 
of deposition rate at the solid surface. 

Thus a new homogeneous Dirichlet boundary condition at the deposition 
surface has been applied to the product species. This boundary condition is 
consistent with the method of calculation of the deposition rate which directly 
uses the diffusion mass flux of the product species at the susceptor surface. In 
the present study product species has been specified as zero at the deposition 
surface with the view that prediction of diffusion flux at the susceptor surface 
will serve as an upper limit to the actual deposition. The reason for this is that 
when a non-negative scalar variable is made to zero at a boundary, the gradient 
of the variable will be maximum. This can be seen in analogy to convective heat 
transfer, where the heat transfer at a surface is maximum when the surface is at 
its coldest possible temperature. 

However, the boundary condition at the substrate surface for all the other 
(reactants) species is just the homogeneous Neumann (no flux) condition. 


3.4 Grid Independence Test 



Figure 3.9: Nondimensional shear stress distribution on the reactor wall for the 
geometry with flat substrate. 
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Figure 3.10: Nondimensional shear stress distribution on the reactor wall for a 
curved substrate. 



Figure 3.11: Nondimensional shear stress distribution on the substrate surface 
for a flat substrate. 

A grid independence test has been carried out in the present study to in- 
vestigate the effect of the grid size on the computed results. The case discussed 
here is the baseline configuration of Re = 100, velocity ratio of the jets V,. = 2, 
angle of the peripheral jet a = 30°, block position Zb = 1- Three grid levels 
namely, 73 x 38, 83 x 48 and 103 x 68 have been studied. Comparative results 
are presented in terms of velocity vectors and shear stress distribution on the 
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Cf 


Figure 3.12; Nondimensional shear stress distribution on the reactor wall for a 
curved substrate. 

substrate surface and the reactor wall. 

Figures 3.9 and 3.10 show the shear stress distribution on the reactor wall 
for the flat and concave substrates. Results for the three grids are quite similar 
in nature for the flat substrate on all the three levels. A slight difference can 
be seen for the curved substrate. Perhaps this is due to the nonuniformity of 
the grid lines and the grid spacing in the annular region between the tip of the 
substrate and the reactor wall. For the three grid levels the peak of the shear 
stress distribution is not at the same axial location in Figure 3.10. The difference 
is not due to grid convergence, but arises from the change in the position of the 
substrate with the number of cells employed in the computation. 

Radial shear stress profiles on the substrate surface have been shown in 
Figures 3.11 and 3.12. For the flat substrate, results compare well with each 
other while there are differences for the concave substrate. As grids are refined, 
a greater amount of information is revealed and the stagnation zone attached to 
the substrate surface is well-resolved. This explains the slightly lower estimation 
of the shear stress on a fine grid. 

The mass transfer code has also been tested for grid independence. Three 
levels of grids have been used for testing the deposition rate of ZnS on the sub- 
strate surface. The results shown in Figure 4.13 match closely with each other 
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Figure 3.13: Radial variation of deposition rate on the solid surface for a fiat 
substrate. 

on all grids. 

Figure 3.14 shows the velocity vectors for the three levels of grids. The out- 
ermost streamline of the recirculation zone has been superimposed in the figures. 
These figures show that results from the three levels of grid are quite similar to 
each other. The position and size of the recirculation patterns are also the same. 
Figure 3.15 shows the contours of the consumption rates of zinc and the mass 
fractions of zinc sulphide. Only the significant contours are shown to maintain 
the clarity of the figures. The results are practically identical. Only the region 
rich in zinc sulphide becomes slightly larger for fine grid. 


It was seen that with an increase in the number of finite volumes, the compu- 
tational time increased approximately as the cube of the number of cells employed 
for the calculations. Thus with a grid level as high as 103 x 68, the economy and 
accuracy of the calculations could not be simultaneously maintained. Due to the 
requirement of a severe computational effort, no attempt was made to generate 
results on a finer grid. Considering the accuracy achieved and the computational 
time involved a moderate grid level of 83 x 48 has been adopted as the basis of 
all subsequent calculations. On the 83 x 48 grid, a flow calculation took twenty 
hours of CPU time; a mass transfer calculation took thirty hours of CPU time. 
A Pentium III, 766 MHz and 512 MB RAM PC was used for all the calculations. 
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Figure 3.14: Velocity vectors for the reactor with a flat substrate, (a) 73 x 38 
(b)83 X 48, (c)103 x 68; Configuration Re = 100, K = 5, = 0.8, Z/, = 1, a = 









Chapter 4 


Results and Discussion 


Results for fluid flow and mass transfer in the CVD reactor are presented in 
this Chapter. In the first section results for the flow and concentration fields are 
presented. The flow field is shown by velocity vectors; significant streamlines are 
superimposed to identify the recirculation pattern, specifically its size, strength 
and location. Contours of the steady state consumption rate of zinc and the 
steady state production rate of zinc sulphide have been shown to get an idea of 
the reaction zone. The mass fraction contours of zinc and zinc sulphide have 
also been shown on a per unit mass of the carrier gas basis. In the following 
section, the deposition rate of zinc sulphide along with the shear stress on the 
substrate are discussed. Flat and curved substrate shapes have been individually 
considered. 


4.1 Flow and concentration fields 

4.1.1 Reactor with a flat substrate 

Figure 4.1(a) shows the flow field for the baseline case (Re = 100,11- = 
2,rb = 0.8, Zj, = 1). At the Reynolds number studied, the central jet reaches 
the substrate surface. The peripheral jet is weak and cannot penetrate into the 
central jet which has a higher momentum. The recirculation pattern is weak 
when the peripheral jet is injected at an angle of 30° towards the central jet 
(Figure 4.2a). This pattern however disappears for the parallel jet configuration. 
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Due to the rapid expansion of both the jets, mixing of the two streams starts 
closer to the inflow plane. A thin stagnation zone appears adjacent to the flat 
substrate. This condition is favourable for uniform deposition. The deposition 
rates would however be low. The uniformity is due to the fact that the carrier 
gas has negligible velocity in the immediate vicinity of the substrate that gives 
zinc sulphide particles enough residence time to be deposited on the substrate 
surface. The concentration fields in Figures 4.1 (b-e) and Figures 4.2(b-e) reveal a 
strong influence of the flow pattern. As the reaction is spontaneous zinc spreads 
a small downstream distance from the inlet plane. Hydrogen sulphide that is 
injected through the central jet spreads throughout the domain. Its distribution 
is also strongly controlled by the basic flow pattern of the carrier gas. The zinc 
consumption rate and the zinc sulphide production rate show that the stretch 
of the reaction zone does not reach the substrate and is confined near the inlet 
plane where the two streams mix. 

Figures 4.3 and 4.4 show the results with the substrate placed at a dimen- 
sionless distance of 0.6 from the inflow plane. The flow fields in the two figures 
are similar, with the exception of an appearance of a recirculation pattern for 
the jet angle a. = 30°. The spread of the reaction zone for the parallel jet case 
(a = 0) is slightly higher due to the delay in mixing of the two streams. The 
zinc consumption zone and zinc sulphide production zone are stretched slightly 
for the parallel jet configuration (a = 0). Zinc spreads a greater distance for 
the parallel jet configuration as in this case the expansion of the peripheral jet 
occurs freely without any constraint. But as the peripheral jet is injected with 
an inward angle it does not expand much and reacts with hydrogen sulphide in 
the central jet closer to the inflow plane. The distribution of zinc sulphide for the 
two jet angles are quite similar. But the regions where maximum concentration 
of zinc sulphide occurs are diff^erent. For the parallel jet configuration, this region 
does not spread to the reactor surface while for the inclined jet configuration, 
a significant amount of ZnS reaches the reactor side wall. The presence of the 
recirculation zone does not have any influence on the concentration distribution. 

Figures 4.5 and 4.6 show the flow field and species distribution when the 
diameter of the substrate is higher, (The blockage ratio being 0.45). A weak 
recirculatory flow occurs for the configuration a = 30° and is completely absent 
for the parallel configuration. It is to noted that with higher blockage levels of 
velocity near the reactor side wall increases and is more uniform compared to 
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the lower blockage geometry. The zinc consumption rate and zinc sulphide pro- 
duction rate clearly show that due to free expansion of the peripheral jet, the 
mixing region is longer for the parallel jet configuration. Hence, the extent of the 
reaction zone is longer for the parallel jet geometry. The distribution of zinc and 
zinc sulphide do not show any significant variation from the baseline calculation 
(Figures 4.1 and 4.2). With higher blockage, a lesser amount of species escapes 
the reaction zone and is indicative of a reduced loss. The recirculatorv pattern is 
once again quite weak to infiuence the species distribution. Hence it has no effect 
on the establishment and size of the reaction zone. 

Figures 4.7 and 4.8 present results for the central-to-peripheral jet velocity 
ratio of 5. Both o; = 0 and 30° cases show strong recirculation. For parallel jets, 
the eye of the vortex appears approximately midway between the inlet and the 
substrate. For a — 30°, the recirculation zone shifts toward the inflow plane. Its 
size being larger, the top of the peripheral jet is reached. The occurrence of a 
strong recirculatory flow can be attributed to the fact that the peripheral jet is 
weak and the main flow envelopes the region above the peripheral jet. As the 
peripheral jet is weaker the extent of the reaction zone is shorter compared to 
the baseline geometry. Thus it can be concluded that the jet velocity ratio has a 
stronger influence on mixing over the jet angle. This is clear from Figures 4.7(b,c) 
and 4.8(b,c). The reaction zones are similar for the two jet angles. The distri- 
butions of zinc and zinc sulphide lead to the same conclusion. Thus for a lower 
jet velocity ratio results are almost independent of the jet angle. A controlling 
infiuence is seen only for higher velocity ratios. 

With a changes in Reynolds number, the flow field and the distribution 
of species are significantly altered. Reynolds number is directly related to the 
mass flow rate of central jet. Figures 4.9 and 4.10 show the results for Re=10. 
The central jet expands in the cross-stream direction faster than in the stream- 
wise direction due to reduced momentum. This improves the mixing of the two 
streams. The size of the stagnation zone adjacent to the flat substrate increases 
as the momentum of the central jet diffuses quickly due to a strong viscous effect. 
Only a weak recirculatory pattern appears for a jet angle of 30°. The reaction 
zones are similar for the two jet angles. This zone spreads a smaller distance 
compared to the baseline case indicating good mixing closer to the inflow plane. 
Distribution of zinc an zinc sulphide show that the species do not spread much 
and accumulate near the inflow plane. The species that is carried by the central 
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jet becomes immobile beyond the central plane of the reactor. Thus zinc sul- 
phide remains in vertical layers in the region adjacent to the substrate surface. 
The species concentration varies near the inlet plane and is nearly uniform at the 
mid-section of the reactor. The effect of the orientation of the peripheral jet is 
suppressed by the strong viscous effect at low Reynolds number. Thus, lowering 
the Reynolds number does not offer any advantage in the reactor performance. 
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Figure 4.1; Velocity vectors (a), steady state contours of consumption rate of Zn 

(b), steady state contours of production rate of ZnS (c), contours of mass 

of Zn (d), contours of mass fraction of ZnS (e), Baseline configuration; Re=100, 

K = 2, r. = 0.8, a = 0,Zb = l 









Figure 4.2; Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, K -%n- 0-8, a - 
inwards, Zh = l 
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Figure 4.3: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re^lOO, V;. = 2, r6 = 0.8, a = 0, 
.6 
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Figure 4.4: Velocity vectors (a), steady state contours of consumption rate ot Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, K = 2, n = 0.8, o; = 30® 
inwards, Zb = 0.6 









Figure 4 5- Velocity vectors (a), steady state contours of consumption rate of Zn 
(b) , steady state contours of production rate of ZnS (c) , contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re-100, V, _ 2, ri, - 0.9, a - 0, 





Figure 4.6; Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), a = 30" 

of Zn (d), contours of mass fraction of ZnS (e), Re-100, r , fe 

inwards, Zb = 'i- 







Figure 4.7: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, Vr = 5, = 0.8, a = 0, 




Figure 4.8; Velocity vectors ( 3 -)} steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, K = 5, rt = 0.8, a = 30 
inwards, Zb = I 




gure 4,9: Velocity vectors (a), steady state contours of consumption rate ^ ^ ^ 

K steady state contours of production rate of ZnS (c) contours of frac«cn 
Zu (d), contours of mass fraction of ZnS (e), Re-10, 1 r 2, n 







Figure 4.10: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=10, Vj- = 2, rb = 0.8, a = 30° 
inwards, Zb = l 





4.1.2 Reactor with a concave substrate 


Results for the CVD reactor with a concave substrate are discussed below along 
with the figures for flow and concentration fields. 

Figures 4.11 and 4.12 show the results for the baseline case of Re = 100, i; = 
2, n = 0.8, Zb = 1. Flow fields for the cases of a = 0 and 30° are practically simi- 
lar. For the geometry with a = 30° a recirculation pattern is visible while for the 
parallel jets case it disappears. The stagnation zone attached to the substrate 
surface is larger when compared to the flat substrate geometry. The recirculatory 
zone occurs just above the peripheral jet and velocity near the reactor wall is small 
everywhere. An exception occurs in the annular region between the substrate and 
the reactor wall. In this region the fluid stream passes through a narrow annulus 
Hence flow accelerates and the velocity is quite high. In the parallel jet geometry 
the peripheral jet freely expands and thus mixing with the central jet occurs over 
a longer distance compared to a = 30°. The reaction zone is not affected due to 
the curvature of the substrate as it is placed at a downstream distance that is 
not enough to alter remarkably the flow field. Thus the stretch of the reaction 
zone is similar to the case of flat substrate. 

The distribution of zinc in the reactor is similar for the two jet angles It 
spreads however for a longer distance for the parallel jot geometry. Zinc sulphide 
spreads over the entire domain up to the substrate and its distribution is iden- 
tical for the two angles. In the present geometry, the concavity of the substrate 
surface alters the flow direction near the obstacle and is not similar to the case of 
a fiat substrate. Flow moves backward near the surface which pushes the region 
rich in zinc towards the inflow plane. Thus the zinc-rich region is smaller and 
nearer to the inlet plane for the concave geometry. This result is not seen in the 
hydrogen sulphide distribution because the central jet through which this species 
is injected has significant inertia in the stream-wise direction. 

Figures 4.13 and 4.14 present the results for the case with substrate placed 
at a dimensionless distance of 0.6 from the inflow plane. There is no significant 
recirculation pattern for the parallel jet case while a weak pattern appears for 
a = 30°. This pattern occurs at an angle to the flow direction for the curved 
block, in comparison to the flat substrate. The size of the stagnation region 
decreases when the substrate is brought closer to the inflow plane. The length 
of the reaction zone is greater for the parallel jet geometry (o: = 0) compaied 
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to a' = 30 . It can be seen in Figures 4.13(b,c) and 4.11(b,c) that the reaction 
zone reaches the reactor side wall to a greater extent when compared to the flat 
substrate, mainly due to its dependence on the basic flow pattern of the carrier 
gas. Zinc sulphide distribution for the curved block also shows a similar devia- 
tion from the flat substrate. The reason is the backflow of the bulk fluid ahead of 
the stagnation zone. Thus for the concave substrate, though more zinc sulphide 
reaches the substrate surface due to its concavity, the reaction zone and the re- 
gion rich in zinc sulphide shift in the upstream direction. 

With the increase in the diameter of the substrate. Figures 4.15 and 4.16, 
the flow field does not change significantly. Trends observable for the flat sub- 
strate geometry are once again reproduced. For the parallel and inclined jets 
ct = 0 and a = 30°, the recirculation pattern is absent. The velocities through- 
out the domain are slightly higher due to a higher blockage ratio. The effect of 
blockage and the concavity of the substrate surface does not lead to a velocity 
distribution that varies significantly in the middle section of the reactor. This 
result was also obtained in the flat substrate simulation. For the parallel jet 
geometry the reaction region is slightly stretched, indicating that mixing occurs 
over a longer distance compared to a = 30°. This reaction zone is quite similar to 
the baseline case. The distribution of zinc is similar for the two geometries with 
different peripheral jet angles. The zone rich in zinc shifts upstream compared to 
the flat substrate geometry. Zinc sulphide distribution also shows a similar trend. 
These observations lead to the conclusion that the blockage ratio has a marginal 
effect on the flow and species concentration fields. 

Figures 4.17 and 4.18 present numerical results when the central-to-peripheral 
jet velocity ratio is 5. For both cases (o; = 0 and a = 30°) a strong recirculation 
pattern occurs near the reactor wall. The size of this pattern however is smaller 
when compared to the flat substrate geometry. This is due to the concavity of 
the substrate which leads to a considerable negative axial velocity. The zinc 
consumption rate and zinc sulphide production rate decrease significantly due 
to incomplete mixing. There are no observable differences in the reaction zone 
between the present case and the one with a flat substrate at a jet velocity ratio 
of 5. Concentration distribution of all the species show that due to lower velocity 
of the peripheral jet, zinc does not spread over an appreciable distance. It leads 
to a reaction zone confined closer to the inlet plane. Moreover, the zone con- 
taining the maximum zinc sulphide shifts to the reactor wall due to the presence 
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of a ^\6ak peiipheral jet. This is certainly an adverse effect with respect to the 
deposition rate and the throughput of the reactor. 

Figures 4.19 and 4.20 show the flow and concentration fields for Re=10. 
Due to larger viscous forces both the central and the peripheral jets expand near 
the inlet plane leading to an early mixing of the streams. The lower levels of 
velocity through out the domain compounded with the effect of concavity of the 
substrate surface results in a large stagnation region. This result was also seen 
for the flat substrate. The absence of convective effects leads to stagnant vertical 
layers of ZnS over a larger portion of the reactor. This distribution is not helpful 
for higher rates of deposition but is favorable if the goal is to form a uniform 
layer. The distribution of zinc shows that due to a higher transverse diffusion 
of momentum this species has been restricted to a limited zone near the inlet 
plane. The region where zinc is consumed and zinc sulphide is produced are also 
restricted to the vicinity of the inlet plane. 



Figure 4.11: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b) , steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Baseline configuration: Re=100, 
Tk = 0 . 8 . a = 0, Zh = ^ 





Figure 4 12- Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, — 2, rj — 0.8, a — 30 

inwards. Zh = l 





Figure 4.13; Velocity vectors (a), steady state contours of consunaption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, = 2, = 0.8, a = 0, 
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Figure 4.14; Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, V = 2, rj, = 0.8, a — 30 
inwards, 2’;, = 0.6 







Kure 4 15: Velocity vectors (a), steady state contours of consumption rate of Zn 
), steady state contours of production rate of ZnS (c), contours of mass fraction 
Zn (d), contours of mass fraction of ZnS (e), Re-100, K - 2, 7 6 - 0.9, o; 0, 





Figure 4.16: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, K = 2, n = 0.9, a = 30'’ 
inwards, Zb = I 


Figure 4.17: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, Vr — 5, rt — 0.8, o: — 0, 










Figure 4.19: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=10, K = 2, = 0.8, a = 0, 
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Figure 4 . 2 O 1 Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=10, == 2, r 5 — 0.8, a — 30 
inwards, = 1 
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4.2 Wall shear stress and deposition rate of zinc 
sulphide 


In this section results are presented in terms of wall shear stress and the 
deposition rate of zinc sulphide on the substrate surface, the dimensionless local 
wall shear shear stress, namely the skin friction coefficient is defined as 


bcaU! \pcaUi 

The deposition rate signifies the volumetric rate at which a zinc sulphide layer is 
added to the substrate. Instead of the absolute thickness, the deposition rate is a 
moie lelevant paiameter since it permits a comparative assessment of the role of 
the flow parameters and reactor geometry. The local deposition rate is calculated 
as follows. 

— = 

dt PZnsAs PZnsA-s 

where s is the thickness of zinc sulphide deposited, is the diffusive mass flux of 
zinc sulphide on the substrate given in Kg/s, is the surface area vector of the 
substrate and DznS is the binary diffusion coefficient of zinc sulphide into argon. 


4.2.1 Properties of a flat substrate 

Figure 4.21 shows the axial variation of the skin friction coefficient along 
the reactor wall, shear stresses are significantly higher on the wail adjacent to the 
annular region. As the bulk fluid stream passes through this narrow annulus fluid 
stream accelerates leading to increase in shear stress. Elsewhere shear stress is 
uniform and quite low in magnitude. It can be seen that values of the skin friction 
coefficient are relatively higher at all the points on the reactor wall for Re=10. 
Shear stress is minimum for jet ratio 5 which shows that at lower peripheral jet 
velocities very small amount of fluid reaches the reactor wall. The variation of 
skin friction coefficient is practically identical for a = 0 and o; = 30° as can be 
seen from the Figures 4.21(a) and (b). Thus effect of jet angles is not felt at the 
reactor wall. 
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Figure 4.22 shows the effect of Reynolds number on the skin friction coeffi- 
cient Cf and the deposition rate of zinc sulphide on the substrate surface. Shear 
stress increases linearly along the substrate surface as one moves away from the 
axis of the leactor. Fairly large shear stresses are encountered at the tip of the 
substrate. Fluid velocities are small in the stagnation zone adjacent to the sub- 
strate. The extent of this region decreases as one progresses towards the tip of 
the substrate, leading to significant velocity gradients. It can be seen that shear 
stress levels are low^er for Re=100. This shows that viscous forces acting on the 
substrate are relatively higher at Re=10. 

The deposition rate is almost invariant closer to the axis and varies signifi- 
cantly around the tip of the substrate. This trend near the edges of the substrate 
is an indication of a nonuniform thickness of the deposited material. The growth 
rate is higher at Re=100 at the tip of the substrate. This is because at higher 
velocities near the substrate surface, particles of zinc sulphide do not have the 
residence time for deposition, leading to a lower deposition rate. In the region 
near the reactor axis deposition rate is small because generation rates of zinc sul- 
phide are negligible (Section 4.1). Results show that jet angle at the inlet plane 
does not have an impact on the deposition rate and the shear stress on the solid 
surface. 

Figure 4.23 shows the effect of the position of the substrate on the shear 
stress and deposition profile on the substrate surface. When substrate is placed 
closer to the inlet plane the fluid stream injected through the central jet reaches 
the substrate leading to higher velocities near the substrate. This leads to a higher 
skin friction coefficient for the case = 0.6 on the substrate surface. This effect 
is significant away from the tip as well. The deposition rate distribution follows 
the trends of shear stress on the surface. The deposition rate at the tip of the 
substrate does not depend strongly on substrate position. The effect of the jet 
injection angle a is insignificant, though the parallel jet configuration predicts a 
marginally higher deposition rate. 

The effect of an increase in the diameter of the substrate is depicted in 
Figure 4.24. Shear stress and deposition rate profiles on the substrate surface aie 
practically identical for the two blockage ratios. With a higher blockage, shear 
stress at the tip of the substrate is approximately twice the value for the lower 
blockage. This is due to the fact that increase in blockage reduces the area of the 
annulus between the substrate and the reactor wall. Thus fluid acceleration is 
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higher compared to the case of lower blockage ratio leading to significant increase 
in velocity gradient in the annular region. Figure 4.21 shows the axial variation 
of the skin friction coefficient on the reactor wall for the increased blockage ra- 
tio. The shear stress is significantly higher at the region above the substrate 
Elsewhere it is almost a constant with a small value. The injection angle of the 
peripheral jet has only a marginal effect on the wall properties. 

Figure 4.25 shows the effect of changing the central-to-peripheral jet veloc- 
ity ratio on the shear stress and growth profile on the substrate surface. The 
interdependence of the distribution of shear stress and the deposition rate on the 
substrate surface for this case is different from the other parameters. A higher 
peripheral jet velocity ( a lower velocity ratio Vj. = 2 compared to 5) shows higher 
levels of shear stress and deposition rate on the surface. This is in contrast to 
the earlier results. The transport aspects of the species are seen to be more im- 
portant than the base flow field of argon. Due to a higher peripheral jet velocity, 
the extent of the reaction zone is larger leading to the production of zinc sulphide 
closer to the substrate. 

For lower velocity ratio, the species front moves very slowly and reactions 
occur within a restricted zone close to the injection plane. Thus the distribution 
of the shear stress on the substrate surface does not give an estimate of the de- 
position profile. The axial variation of the shear stress on the reactor wall follows 
the same trend over a range of parameters. When the velocity ratio is higher, 
a slightly higher shear stress in the tip region of the substrate is observed. The 
injection angle of the peripheral jet does not have any significant effect on the 
reactor properties. Thus changing the jet velocity ratio from 2 to 5 at Re=100 
does not improve the deposition rate quantitatively but becomes more uniform 

along the substrate. 
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Figure 4.21: Axial variation of the skin friction coefficient on the wall of the 
reactor tube for the fiat substrate geometry; ck = 0 (a), a = 30° (b). 
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Figure 4.22. Effect of Reynolds number on the skin friction coefficient and the 
deposition rate on the substrate surface, n = 0.8, V; = 2, = 1. 
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Figure 4.23: Effect of the position of the substrate on the skin friction coefficient 
and the deposition rate on the substrate surface, Re=l00, = 0.8, VJ. = 2. 
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Figure 4.24: Effect of the diameter of the substrate on the skin friction coefficient 
and the deposition rate on the substrate surface, Re=100, K = 2, Zj, = 1. 




Figure 4.25: Effect of the central-to-peripheral jet velocity ratio on the skin 
friction coefficient and the deposition rate on the substrate surface, Re=100, 
n = 0.8, Zb = 1 . 
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Figure 4.24: Effect of the diameter of the substrate on the skin friction coefficient 
and the deposition rate on the substrate surface, Re=100, Vr = 2, Zi, = 1. 



Figure 4.25: Effect of the central-to-peripheral jet velocity ratio on the skin 
friction coefficient and the deposition rate on the substrate surface, Re=100, 
n = 0.8, Zb = 1 . 
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4.2.2 Properties of a concave substrate 


The shear stress distribution on the solid surfaces and the deposition char- 
acteristics on a concave substrate have been presented below. In the figures, the 
effect of various parameters has been depicted by varying one of them at a time 
with respect to the baseline case (Re = 100, K = 2,rb = 0.8, Zb = l). 

Figure 4.26 shows the axial variation of skin friction coefficient on the re- 
actor wall. Levels of shear stress are low and uniform along the reactor wall 
A sudden irrerease in shear stress occurs above the substrate surface. This re- 
gion moves upstream when compared to the fiat substrate configuration. This 
is because the annular gap where fluid accelerates is located approximately at a 
dimensionless axial location of 0.4 from the inlet plane. 

Reynolds number does not drastically alter the skin friction coefficient and 
the deposition rate on the substrate surface as seen in Figure 4.27. Although at a 
lower Reynolds number, viscous effects predominate over inertia forces, the skin 
friction coefficient does not show any significant variation along the substrate 
surface with Reynolds number. It is only near the tip the substrate surface that 
the viscous effect is felt and shear stresses are relatively higher for Re=10. 

There is a reduction in the shear stress at the tip of the substrate at both 
Reynolds numbers. The result does not match the flat substrate profile. One 
reason for this trend could be the numerical treatment of the corner cells. As 
zero penetration of mass has been implemented as the boundary condition, local 
tangential velocity is nonzero at the cell faces. This leads to a reduction in the 
velocity gradient and a lowering of the shear stress. 

Figures 4.28 shows the distribution of skin friction coefficient and the de- 
position rate on the substrate surface for a block position of Zb = 0.6. As the 
substrate is placed closer to the inlet plane, the velocity front from the central jet 
completely reaches the substrate and the size of the stagnation zone decreases. 
This leads to higher velocity near the substrate and consequently higher shear 
stresses. Due to the concavity of the substrate surface, appreciable negative axial 
velocity occurs when the bulk of the fluid is diverted along the substrate. The 
evolution of the negative axial velocity retards zinc sulphide particles in their 
movement towards the substrate surface. The stagnation zone in the immediate 
vicinity of the surface helps the diffusing species to overcome the convective effect 
and deposit on the substrate. It should be noted that the growth rate profile and 
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the shear stress distribution at the tip of the substrate follow closely the trends 
for a flat substrate. This emphasizes the role of species diffusion in deposition 
with respect to the bulk convection of the fluid. 

Figure 4.29 shows the effect of increasing the diameter of the substrate and 
hence the blockage ratio. Both skin friction coefficient and the deposition rate on 
the substiate surface do not vary significantly along the substrate surface. This 
is siniilai to v hat happens in the case of a flat substrate. Only a slight deviation 
of the sheai stiess and growth rate is observable at the tip of the substrate. This 
lesult has also been discussed in context of the flat substrate. A minor variation 
in the shear stresses and the deposition profile along the substrate surface for 
different blockage ratios shows that increasing the diameter of the substrate does 
not alter the flow field or the deposition characteristics. Figure 4.26 shows the 
axial variation of the skin friction coefficient on the reactor wall for a blockage 
ratio. At a higher blockage, the shear stress is higher due to an increase in the 
velocity gradient at the substrate tip. This occurs due to the acceleration of the 
fluid stream passing through the narrow annular gap. 

The effect of the central-to-peripheral jet velocity ratio on the shear stress 
distribution and the growth rate are depicted in Figure 4.30. No deviation in the 
shear stress and the deposition rate is observed with respect to the flat substrate 
configuration for a given set of parameters. The deposition profile for the two 
velocity ratios are identical in nature. Due to the longer reaction zone at a higher 
velocity ratio, a greater amount of zinc sulphide reaches the substrate, thus in- 
creasing the average level of deposition. The axial distribution of the skin friction 
coefficient on the reactor wall for a higher velocity ratio is shown in Figure 4.26. 
The values of shear stress are nearly equal everywhere except in the annular gap 
at the substrate location. With a higher peripheral jet velocity a greater amount 
of flow occurs through the annulus leading to a slight increase in the velocity 
gradient. 
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Figure 4.27: Effect of Reynolds number on the skin friction coefficient and the 
deposition rate on the substrate surface, n = 0.8, V; = 2, Zt, = 1. 



Figure 4.28: Effect of the position of the substrate on the skin friction coefficient 
and the deposition rate on the substrate surface, Re=100, rt, = 0.8, = 2. 






radial coordinate radial coordinate 




■p 0.25 


0 0.05 0.1 0.15 0 2 0.25 0.3 0.35 0.4 0.45 
Cf 


rb=0.9, a=30° -A 
rj 3 = 0 - 8 , a=30° -o- 
rb=0.9, a=0 -X- 
rj 5 = 0 . 8 , oc =0 — 


deposition rate (|a.m/min) 


Figure 4.29: Effect of the diameter of the substrate on the skin friction coefficient 
and the deposition rate on the substrate surface, Re=100, K = 2,Zb = 1. 
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Figure 4.30: Effect of the central-to-peripheral jet velocity ratio on the skin 
friction coeflBcient and the deposition rate on the substrate surface, Re— 100, 
n = 0.8, Zh = 1 . 






4.3 Evaluation of alternative geometries 

lu U.is scctio.. limited results for the alternative geometries have been pre- 
sentec . T^re schomatm views of these geometries are shown Figures 1.6(a) and 
1.0 b) of Chapter 1. The first of these geometric consists of a substrate whose 
suifaM rs comnx facing the main flow direction. The second has coaxial jets 
with H..S and Zi, are injected through the central and the annular nozzle respec- 
tively. It IS expected that when the convex surface of the substrate faces the 
flow significant changes m the flow distribution occur, in turn leading to different 
deposition characteristics. The coaxial geometry is expected to extend the size of 
the reaction zone towards the substrate and once again improves the deposition 
characteristics of the CVD reactor. 


4.3.1 Flow and concentration fields with a convex sub- 
strate 

In this section the flow and the concentration fields hnve been discussed 
with the parameters set at Re = 100, K- =^2, n = 0.8, = 1 and a = 30°. The 

significant contours of species distribution, their consumption rates and produc- 
tion rates have been shown. This approach has retained the clarity of the figures. 

The baseline case employed for comparison of different parameters is iden- 
tical to that of the concave substrate configuration (Re = 100, Vr = 2, = 

0.8, Zb = 1 and a = 30°). Figure 4.31 shows the flow and concentration fields 
for the baseline configuration. Due to the convexity of the substrate surface the 
stagnation zone is seen to disappear. The geometric shape of the substrate aids 
the fluid stream to reach closer to the solid surface. This is an important dif- 
ference with respect to the previous geometries. The expansion of the central 
jet occurs up to a greater distance downstream compared to the previous geome- 
tries. Moreover, uniformity of the velocity at the central zone of the reactor is 
improved, leading to enhanced transport of the species. The shape and size of 
the recirculation pattern however do not change. This is because the presence of 
the peripheral jet that confines it in the cross-stream direction, restricting it to 
the corner of the reactor. The flow field near the inlet plane remains unaltered as 
the geometry of the substrate does not propagate sufficiently upstream against 
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the inomcntum of the central jet. This leads to a reaction zone of almost the 
same length as the reactor with a concave substrate. The distribution of zinc 
is similar to that of the concave block due to the absence of any changes in the 
flow field near the inflow plane. The region rich in zinc sulphide is located at the 
same position and its size and shape also do not change. This emphasizes that 
at low Reynolds numbers it is the basic flow field that controls the transport of 
the species. Changes in the shape of the substrate do not alter the flow field due 
to the presence of the strong central jet. 

Figures 4.32 shows the flow field and species distribution for the case of 
substrate being placed closer to the inflow plane. The stagnation zone reduces 
even further compared to the case of Z;, = 1. Expansion of the central jet occurs 
quite smoothly due to the absence of negative axial velocity near the substrate 
surface. Levels of velocity near the substrate are higher compared to the two ear- 
lier geometries. Severe non-uniformities in the velocity profile in the central zone 
observed in the case of a concave substrate are absent here. This is due to the 
convexity of the substrate that allows the fluid stream to pass by the substrate 
without moving upstream. The recirculation zone does not change in shape and 
size with respect to the block position of = 1. Upward stretching of the re- 
action envelope towards the annular gap here is a favorable feature. Thus the 
transport of the species occurs in the stream-wise direction towards the substrate. 
The distribution of zinc shows that the species-rich regions do not get compressed 
by the backflow of the bulk fluid. This adverse effect was weakness of the con- 
cave substrate. The zinc sulphide-rich zone increases in size and is a continuous 
layer because the reduced upstream diffusion of momentum does not interrupt 
the mixing. This shows smooth mixing of the two fluid streams injected through 
the central and peripheral jets for the present geometry has been accomplished. 

Figure 4.33 shows the results for a central-to-peripheral jet velocity ratio 
of K = 1 Velocity levels inside the reactor increases while the cross-stream ve- 
locity is small everywhere. This is due to the fad that the fluid stream injected 
through the peripheral jet has a higher momentum and flows in the axial direc- 
tion after meeting the central jet. It can be noted that the mixing interface is 
fiat and spreads in the stream-wise direction. This region also advances towards 
the substrate leading to an enhanced transport of the species. With the increase 
in the velocity of the peripheral jet. the reaction envelope stretches to almost 
twice the distance compared to the higher jet velocity ratio. The distribution of 
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zinc sulphide is also favourable for deposition compared to a higher peripheral jet 
velocity. Due to the longer reaction zone in the stream-wise direction, the shape 
of the zinc sulphide-rich zone also changes. This zone is increasingly parallel to 
the react.or axis compared to the case of a lower peripheral jet velocity where it 
spreads in the cross-stream direction. This is a feature which is desirable since 
the growth rate is expected to increase and results in lower material losses. 
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Figure 4.32: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, K = 2, r(, = 0.8, a = 30° 
inwards, Zh = 0.6 









Figure 4.33: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b) , steady state contours of production rate of ZnS (c) , contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100, V = !> ^6 = 0-8, ck = 30 
inwards, Zb = l 






4.3.2 Flow and concentration fields with coaxial jets 


Foi the coaxial jet configuration the parameters considered are: Re — 
100, Iv = 1,2, 7'i = 0.8, Zh = 1,0.6. Calculations have been carried out with a 
conca,vc substrate. 

Figure 4.34 shows the flow and species concentration fields for the baseline 
case of Re = 100, = 2, = 0.8, = 1. The central jet does not expand 

freely due to the confining effect of the peripheral jet. Mixing starts at the inlet 
plane as the jets meet after their injection. The recirculation region is stronger 
as the peripheral jet has been moved towards the axis of the reactor. This region 
is at a small angle with the axial direction, thus revealing a gradual expansion 
of the peripheral jet. The strong recirculating flow brings a certain portion of 
the zinc stream to the vicinity of the inlet plane. But as the initial zinc stream 
envelopes the central jet only a fraction of the recirculated zinc particles react 
wdth hydrogen sulphide. The shear between the two jets at their interface leads 
to vigorous and accelerated mixing. The reaction zone is thus limited to a small 
distance from the inlet plane. A secondary zinc sulphide production zone is ob- 
served at the boundary of the peripheral jet and is attached to the inflow plane. 
This is due to the reaction between zinc particles moving in the recirculation zone 
and traces of hydrogen sulphide available in the vicinity of the inlet plane. Zinc 
is distributed in the region very close to the inlet plane, the reason being th® 
accelerated mixing and spontaneous reaction of the injected species. The regions 
of large zinc sulphide which show the zone of mixing with the most favourable 
stoichiometry is almost reduced to a point. This again emphasizes the adverse 
role of the accelerated mixing and reaction on the zinc sulphide distribution. 

Figure 4.35 shows the flow and concentration fields for the case of coaxial 
jets when the substrate is placed closer to the inflow plane. The shape of ti 
recirculation pattern is now different and is closer to being circular. The 
of concavity of the substrate surface is felt by the fluid stream at a plane 
to the injection surface and the expansion of the peripheral jet is impeded 
inflexion point in the velocity profile is clearly visible near the mid plane o t re 
reactor, x/d = 0.3. The shearing effect of the two streams at their inter ace 
occurs over a shorter distance due to a closer placement of the substrate rom 
the inflow plane. The reaction envelope is similar to the baseline case, where 
Zft = 1. The secondary zinc sulphide production zone however increases i 
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aiul mixes with the primary zone. Also a secondary zinc consumption zone is 
vdsiblc. This shows that the recirculatory flow is stronger in the present case 
owing to a sharper variation of the axial velocity in the cross-stream direction 
compared to the baseline case {Zb = 1). The distribution of zinc shows that 
the zinorich zone is compressed owing to greater amount of reaction with hy- 
drogen sulphide. The maximum zinc sulphide-rich zone stretches towards the 
reactor axis. These zones were practically a point for the baseline case. Thus 
wdth the substrate placed nearer to the inlet plane, strong recirculating current 
increases the reaction between the trace quantities of hydrogen sulphide and the 
recirculated zinc in the secondary mode, leading to a larger volume where the 
zinc sulphide concentration is quite high. 

Figure 4. 36 shows results when the central-to-periphcral jet velocity ratio 
Vr is reduced to unity. The recirculation pattern is similar in shape, size and 
position to the baseline case of K = 2. The shearing effect at the jet interface 
decreases as their velocities are now equal. A velocity profile at any cross-section 
is smoother, and strong gradients are absent. The reduction in shear between the 
two streams leads to gradual mixing followed by reactions This is manifested in 
the increased extent of the reaction zone. Zinc consumption rate and m turn zinc 
sulphide production rate significantly increases. Quantitatively they are almost 
five times the rates for the baseline case of 14 = 2. The gradual mixing of the 
two streams leads to a reaction zone that spreads in the stream-wise direction 
pointing towards the substrate. The region rich in zinc increases due to the higher 
peripheral jet velocity as well as its free expansion. This is aided by the absence 
of a strong backflow that occurs at a block position of Zb — 0.6. The size 
maximum zinc sulphide region also increases due to the gradual mixing of the 
streams over a longer downstream distance from the inflow plane. 
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Figure 4.34: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c), contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Baseline configuration: Re-100, 

Vr = 2. Th = 0.8, a = 0, Zb = 1 
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Figure 4.35: Velocity vectors (a), steady state contours of consumption rate of Zn 
(b), steady state contours of production rate of ZnS (c) , contours of mass fraction 
of Zn (d), contours of mass fraction of ZnS (e), Re=100j = 2, = 0.8, a = 0, 
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4.3.3 Wall shear stress and Deposition rate Profiles for a 
convex substrate 


Distributions of the skin friction coefficient and the deposition rate on the 
substrate surface have been presented in this section. Also skin friction coefficient 
on the reactor wall has been shown. 

Figure 4.37 shows the axial variation of the skin friction coefficient on the 
reactor w'all. The trends are same to those of earlier geometries (Figures 4.21b 
and 4.26b). The values of Cj is considerably lower for the present case. The 
reason for this is the uniformity in the flow field. Large velocity gradients that 
was obtained for the fiat and concave substrate geometries are absent here. The 
convex shape of the substrate has a streamlining effect on the flow which passes 
over it without moving upstream. 

Figure 4.39 shows the radial variation of the skin friction coefficient and 
the deposition rate on the convex substrate for three sets of parameters. Shear 
stresses on the substrate surface are comparable for the three cases. Near the 
axis the shear stress is slightly higher due to higher velocity levels near the sub- 
strate surface when the block position is Zb = 0.6. It should be noted that shear 
stress increases at the tip of the substrate, a result that is the reverse of what 
was presented for the concave surface. The explanation for this anomaly is the 
following. For the convex substrate the stagnation zone completely disappears 
leading to a considerably large velocity in the vicinity of the substrate. Thus the 
velocity gradient is also greater in magnitude along the surface. At the tip of 
the surface, the fluid stream transits from a gradual passage to a narrow annulus. 
This leads to acceleration in the fluid stream. It is to be noted that in the present 
configuration the tip of the substrate lies near the annulus. This is not the case 
for the concave substrate due to its finite thickness. Thus the proximitj o 
substrate tip and the narrow annulus results in a greater velocity gradient and 
in turn the skin friction coefficient. The effect of fluid acceleration completely 
overwhelms the presence of a low tangential velocity related 
^ling of the corner cells in the numerical algorithm. This issue has ^een ^ ss d 

r the context of an increase of shear stress at the tip o the -c" 

. . . flip trpTid of shoEi stross on biic 

ffie distribution of the deposition rate follows the tre 

1 -* u frirfinTi alons substra.t6 tor a dIock 

ubstrate surface. Due to a higher skin friction along ,„K,trate 

osition of Z. = 0.6, the deposition rate is low at all po. 
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With a ceutral-to-poripheral jet ratio Vj. = 1, the growth rate is considerably 
liigher coniparod to the other case of Vj. = 2. This occurs due to the stretching 
ot the reaction zone in the axial direction as a greater amount of zinc sulphide 
is transport.(*d to t,he substrate. There are no significant deviations in the axial 
distribution of the skin friction coefficient on the reactor wall with reference to 
the concav(' substrate geometry. This shows that the effect of the shape of the 
substratr^ is not, felt in the region near the reactor wall. 


4.3.4 Wall shear stress and Deposition rate Profiles for 
coaxial jets 


Figure 4.38 shows the axial variation of the skin friction coefficient on the 
reactor wall. The trends are similar to the earlier case (Figure 4.37). The peaks 
of the shear stress distribution are now different. Due the changes in the di- 
rection of curvature of the substrate surface, the position of the annular region 
between the substrate tip and the reactor wall is different for the two cases. Thus 
maximum shear stress occurs at a different axial location for the two geometries 
(Figure 4.37 and 4.38). The values of (7/ for the present case is smaller than 
the cases considered in section 4.2. This is due to the fact that for a coaxial jet 
configuration both the central and the peripheral jets constitute the bulk flow, 

leading to a uniformity in the velocity field. 

Figure 4.40 shows the variation of the skin friction coefficient and the depo- 
sition rate on the substrate surface for the coaccial jet configuration. The slope of 
the shear stress along the surface is higher for the three sets of parameters studied. 
This is due to a greater uniformity of the Sow field. Skin friction coefficient is also 
higher for the case when the substrate is placed closer to the inflow plane. In this 
configuration, the streams from the central and the peripheraljet constitute the 
bulk flow. As a greater amount of fluid stream passes over the substrate that has 
a poor streamline shape, an appreciable amount of shear and hence skin friction 
occurs along the surface. The values of the skin friction coefficient are compa- 
rable with the configuration when the jets are separated by a certain drst^ca 
A reduction in the shear stress is observed at the tip of the substrate has 

been pointed out in the discussions for the concave -^atrate “d je ts th^are 
not coaxial. For the present simulation, the deposition rate ^ ^ 
quite uniform (Figure 4.40). But the increase in the growth ra e 
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substr&t6 continues to occur. 'This feature has also been observed for the concave 
substrate with separated jets. The deposition rate of zinc sulphide is higher at 
all the points on the substrate due to the axial alignment of the reaction zone. 
The deposition rates for the other two sets of parameters are also close. Thus 
for a coaxial jets configuration, the uniformity of the central and peripheral jet 
velocities is favourable for a uniform as well as higher levels of the deposition rate 
The axial variation of the skin friction coefficient on the reactor wall is shown in 
Figure 4.38. The trend is similar to the earlier geometries studied in section 4.2. 
Only a slight variation is observed in the region between the inlet plane and the 
substrate. The variation is due to the backflow near the reactor tube wall for the 
present geometry. 
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Figure 4.37; Axial variation of skin friction coefficient on the reactor tube wall 
for the convex substrate geometry; a = 30°. 
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Figure 4.38: Axial variation of skin friction coefficient on the reactor tube wall 
for the coaxial jets configuration; 0 = 0. 
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Figure 4.39: Skin friction coefficient and the deposition rate on the substrate 
surface for the convex substrate geometry with Re=100, a = 30° and n = 0.8; 
(a) K = 1, Zb = 1; (b) Vr = 2, Zb = 0.6; (c) K = 2, Zb = 1. 



Figure 4.40: Skin friction coefficient and the deposition rate on the substrate 
surface for the coaxial jets configuration with Re=100, a = 0 and rj = 0.8; 
(a) K = 1, Zb = 1; (b) K = 2, Zb = 0.6; (c) K = 2, = 1. 






Chapter 5 


Conclusions and Scope for Future 
Work 


Major conclusions arrived at in the present work are summarized in the present 
Chapter. 


5.1 Conclusions 


The fluid flow and mass transfer phenomena inside a low pressure CVD reactor 
have been studied via a mathematical model. Of special interest is the deposition 
characteristics of zinc sulphide on a passive substrate from zinc and hydrogen 
sulphide in the gas phase. The results obtained in the present study help judging 
the merit of a reactor when the geometric configuration and flow parameters are 
changed. The following conclusions have been arrived at in the present work 

1. For a flat substrate placed at right angles to the reactor axis, the flow field 
mimics the properties of a two-dimensional stagnation flow. Here the streamlines 
as well as the bulk of the product species bypass the substrate. Hence, the 
deposition rate on the substrate is uniform but the magnitude of the deposition 
rate is low. 

2. When a concave surface is employed, the stagnation region is further 
enlarged over which the fluid velocities are small. The bypassing of the concave 
substrate by the product species continues £is in the flat substrate. Once again, 
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the d('pc)sitii)ii rate is uniform but is quite small to be of practical value. 

3. A convc'x substrate is better streamlined, and allows the bulk of flow and 
the ivact iii}; species t-o come close to it. This results in a much smaller stagnation 
zone. Thus con\-ecf.iv(' flow significantly contributes towards the deposition of zinc 
sulphide on the sul)strate in addition to species diffusion. Hence the deposition 
rate iiuu'cases considerably over and above what is achieved by concave and flat 
substrates. The average deposition rate for a convex substrate was computed 
to be 1.027 /im/min, compared to 0.542 /rm/min for a concave substrate, an 
improwment of a factor of two. The degree of uniformity of the deposition 
was also found to be superior for a convex geometry. Unfavourable edge effects 
endemic to flat and concave substrates were also mitigated. 

4. If the incoming jets are arranged coaxially, the streams of reactants mix 
and react at their interface. The moving zinc sulphide front reaches the substrate 
surface in a manner that favours a uniform and a higher deposition rate. The 
presence of a large recirculatory flow around the jets brings zinc particles closer 
to the unmixed hydrogen sulphide particles, thus improving the reaction zone. 
An average deposition rate of 0.527 jum/min was observed for the coaxial jet 
arrangement. Benefits seen here are in terms of the uniformity of deposition and 
a simplification of the technology itself. 


5.2 Scope for future work 

The present research has opened up new possibilities for depositing uniform zinc 
sulphide layers at a greater speed. These trends should be explored through a 
full scale simulation. Factors to be examined are: 

1. Effect of higher Reynolds numbers. 

2. Block position from the inflow plane. 

3. Jet inclination 

Growth of other infrared sensor material such as zinc selenide (ZnSe) can also be 
explored using the model developed. 
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Appendix A 


Evaluation of species properties 


Calculations of the transport properties of the individual gases present in 
the reactor as a mixture have been presented below. The treatment is valid 
for diffusion of a particular species in argon, with the presence of other species 
neglected. This is the approximation of binary diffusion, valid at the limit of 
the dilute solution. All species are in the gaseous state and assumed to follow 
the ideal gas law. This approximation may not be strictly valid for zinc sulphide 
(ZnS) which is formed as a finely dispersed solid phase during the reaction. 


A.l Binary diffusion coefficients 


Values of binary diffusion coefficients at low pressure such as 0.03 atm 
and high temperature such as 1000 K are not readily available in the pub- 
lished literature. In the present study these coefficients have been calculated 
by Chapman-Enskog formula that is in turn derived from the kinetic theory of 
gases Hirschfelder et al. (1964) and Bird et al. (1960). The formula is given as: 


Vab = 1-858 10"^ 



PaAB^ 


where Vab is the diffusion coefficient of species A in B in units of cm“-^ sec~\ T 
the local temperature is in Kelvin, P is the bulk pressure in atmospheres, ct^b and 
^v,AB are the Lennard-Jones parameters which depend upon temperature and 


132 



the mter-molecular potential distribution between the two species. The collision 
integral has been tabulated against the values of ^ for different pairs of 

gases by Bird et al. (1960). The values required to calculate the binary diffusion 
coefficients between different pairs of gases has been given in Table A.l for a 
mean reactor pressure P = 0.0296 atm and a temperature T = 917 K, typically 
seen at the mid-level of the reactor. 

When Lennard- Jones force parameters are not available, they can be calcu- 
lated from the critical property data Hirschfelder et al. (1964) as 

K ^ = 0.757; K 

Here W and Tc are the critical volume ( — ^^j ) and temperature (K) of the 
species for which these constants are to be calculated. The values of binary 
diffusion coefficients in the present study have been calculated as: 

^H2S-A = 3.848 10“^ mV\ T^Zn-A = 3.282 10"^ mV^ 

^ZnS-A = 2.457 10"^ m^s"', Pu -A = 1-753 lO'^ m^s"' 

jL 


A. 2 Kinematic viscosity 


Kinematic viscosity of argon given in Pa s has been calculated from the 
following temperature polynomial 

fji = 22.7xl0~®-f7.92xl0~® r+2.93xl0"^^ T^-1.73x 10“^^ T^-1.42x 10"^® 


at the mean reactor temperature, T = 917 K. 


Table A.l: Leonard Jones force parameters for different species 


species 

molecular weight 

a{k) 

f(K) 

A 

39.94 

3.418 

124 

H 2 S 

34.08 

3.733 

221.1 

Zn 

65.38 

2.673 

2377.5 

ZnS 

97.38 

3 

4037.6 

Hs 

2.016 

2.915 

38 
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Appendix B 


Validation of Fluid Flow and 
Mass transfer Codes 

B.l Fluid Flow 


The fluid flow code used in the present study has been validated against the 
results of Seider and Churchill (1971). These authors have studied the mixing 
of conflned jets in the entrance region of a tubular reactor. The equations for 
the conservation of momentum and mass were solved numerically. The physical 
region downstream of the exit of an axial jet mixing with an annular jet of the 
same fluid was considered. In the investigation uniform temperature and uniform 
fluid properties were assumed. Two mixed parabolic-elliptic equations for the 
stream function and vorticity were solved by an implicit finite-difference method. 
Results were obtained for diff’erent jet velocity ratios, jet diameter ratios and 
Reynolds numbers based on the average velocity at the inlet and the diameter 
of the reactor. Figures B.l (b), (c) and (d) show the comparison between the 
computed results of the present work and those reported by Seider and Churchill 
(1971). Axial velocity profiles at three sections have been compared. The two 
sets of results are seen to match very well. 
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B.2 Mass Transfer 


Most of the previous numerical and experimental studies on the transport 
modeling inside a CVD reactor have been directed towards epitaxial growth of 
silicon with surface reactions. Thus benchmark results for mass transfer are not 
available. In the present study the mass transfer code has been validated by 
comparing the results obtained by using the operator splitting algorithm (OS) 
against the results obtained by solving the full transport equations. Time history 
for H 2 S and ZnS mass fraction at three points inside the reactor have been shown. 
Point 1 is located on the reactor wall at an dimensionless axial distance of 0.6 
from the inlet plane. Point 2 is located on the substrate surface at a dimensionless 
radial distance of 0.15 from the axis of the reactor. Point 3 is in the fluid flow 
region located at a radial and axial distances of 0.15 and 0.3 respectively, from 
the inlet plane. 

Figures B.2(a-f) show the transient variation of H 2 S and ZnS at points 1, 
2 and 3 respectively. Results for H 2 S match exactly at all the points for the two 
solution strategies. The transients of ZnS are slightly different at all the points. 
It can be seen that discrepancy in the ZnS solution from the two approaches 
occurs at higher times and is a maximum at steady state. 

Both the numerical techniques utilize first-order accurate time discretization 
and so the numerical error incurred are expected to be comparable. Time step 
were kept equal for the two approaches. As the operator splitting algorithm uses 
a partial analytical solution it is expected that the associated solutions are more 
accurate. Hence, results obtained by the OS algorithm alone have been presented 

in Chapter 4. 
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Figure B.2: Transient variation of species; hydrogen sulphide (a) point 1, (c) 
point 2, (e) point 3; zinc sulphide (b) point 1, (d) point 2, (f) point 3. 










Appendix C 


Derivation of Orthogonal grid 
generation equation 


The equations for generating orthogonal curvilinear coordinates are written 
in Chapter 3. Derivation of theses equations are outlined here. 


C.l Definition 


Definitions of important terms needed in the derivation of the orthogonal 
mapping equation are presented below Thompson et al. (1985). 

A. Covariant base vectors and Covariant metric tensor 


Tangent vectors to the three coordinate lines are the covariant base vectors of the 
curvilinear system. They are designated by 

ai = r^i = 3) (C-l) 

The covariant metric tensor finds its origin in the evaluation of an incremental 
arc length along a spatial curve given by 

{dsf = ai.a.jd^W =gijdee (C-2) 
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The .nere.nenta, arc eeegu. a„ 

terms of the dot products a..aj (j,j = 1,2,3), which forms a symmetric second 

order tensor. These quantities are the components of the covariant metric tensor, 

given by 

ptj = a^.aj {i,j = 1,2,3) 


B. Contravariant base vectors and contravariant metric tensor 


The normal vector to a coordinate surface on which the coordinate f is constant 
is given by 


= (z = l,2,3) 


(C.4) 


The three normal vectors to three coordinate surfaces are the three contravariant 
base vectors of the curvilinear coordinate system. The components of the con- 
travariant metric tensor are the dot products of the contravariant base vectors 

^*^ = aha^ (i,j = 1,2,3) (C.5) 


The relation between the contravariant base vectors and covariant base vectors 
is given by Thompson et al. (1985): 


a' = i(aixa.) 


{i = 1,2,3; ij,k cyclic) 


{C.6) 


with g being the Jacobian of the transformation, namely 

\/? = \/ 9 \ i 9229 zz (C-7) 


C. Derivative Operators 


Expressions for the derivative operators, such as gradient, divergence, curl and 
Laplacian are obtained by applying the divergence theorem on a differential vol- 
ume increment bounded by the coordinate surfaces Thompson et al. (1985). The 
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expression tor aivergence is given in conservaiive lorm a* 


1 

V.A = — ^[(aj X (i, j, k cyclic) 

2=1 


(C^8) 


and in non conservative form, 


V.A = — ^[(a^ X aA;).^f.] (z, j, k cyclic) 

V ^ 2=1 


(C.9) 


with A being a tensor of any order (> 1) 


Using the relation (C.1)-(C.9) we can derive the Laplacian operator as : 

2=1 J = 1 J = l 


C.2 Elliptic generation system 


The simplest elliptic system that satisfies the extremum principle and the 
smoothness property is the Laplace equation given by 

V2f = 0 (i = l,2,3) (C.ll) 

Equation (C.ll) can be applied to each of the coordinate lines and in 

three dimensions. As the coordinate system in the transformed region are made 
of contiguous rectangles its easy to treat the curvilinear coordinates (f) as the 
independent variables and physical coordinates (x^) as the dependent variables in 
the numerical computation. The transformation of equation (C.ll) is obtained 
using equation (C.IO) leading to 

= (C.12) 

i=l j=l j=l 
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where r is the position vector ot a point in me paysiuai uuiudiu, wiiubcui ds 
r = , x"" , x'^) . Since V^r = 0 by definition, we get the elliptic generation 
equation with the help of equation (C.ll) as 

3 3 

(C.13) 

7=1 

C.3 Orthogonal generation system 


The characteristic criterion for orthogonality of the body-fitted coordinates 
is the vanishing of the off-diagonal elements of the covariant and contravariant 
metric tensor, namely 


g,j = ff” = 0 for ! # i 


(C.14) 


For brevity, we write the diagonal elements of the covariant metric tensor as 


K = (! = 1,2,3) 


(C,15) 


From equation (C.12) the Laplace equation for grid generation in terms of the 
physical coordinates as dependent variables is 

1=1 j=i j=i 

On applying orthogonality constraint, equation (C.16) becomes 

= 0 (C.17) 

t=l !=1 
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Ihc first term of equation (C.17) can be written in terms oi covariani metric 
tensor elements as follows, i, j, k being cyclic : 


= a'.a* 

= -^[(aj X X a*)] 

y 

= -[(%-aj)(a/,.afc) - (aj.afc)(a^.afc)] 

= ^[(araj)(aA.ai)] 

y 

1 

~ ~ QjjQkk 

y 

^ QjjQkk 

\/9 y/ 9ii9jj9kk 

__ 1 ^/OjjQkk 

\/9 \f^i 

1 hj /l/j 

y/g hi 

Therefore 

(* = 1. 2, 3)(f, j, k cyclic) 

The second term of equation (C.17) can be shown to be 

[ij^k cyclic) 

y9 hi 


Hence equation (C.17) becomes 


V5tl 






Equation (C.20) can be written in compact form as 


X^(^T^rfOe = 0 (^. 3^ ^ cyclic) 


i=l 


hi 


(C.18) 


(C.19) 


(C.20) 


(C.21) 
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Equation (C.21) is the generating equation for orthogonal mapping, ine two 
dimensional form of this equation used in the present work is written as; 


dC'' dC dr,'- f 8rt' 


dC dC dr fdri’ 


(C.22) 

(C.23) 


wdiere / is the distortion function given by 

J. _ \/g22g33 _ \/^l + 
\/^ + 


{CM) 


with ^33 = 1 for two dimensional geometries. 
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